PCA-HMM (hard regimes, MV per-label, “paper” windows)¶
#!/usr/bin/env python3
# Test 2/5 — PCA-HMM (cross-asset structure, hard regimes) — FAIR VERSION + SAFE TURNOVER PLOT
# -------------------------------------------------------------------------
# Your requests:
# - Keep original (non-causal) regime weights estimation (uses full-sample labels) ✅
# Kept fixes:
# - Constant-vol Risk-Parity baseline (IVP, rolling)
# - Risk-matched (target-vol) headline table @10% vol
# - Net-of-costs with turnover + cost grid (5/10/25/50 bps)
# - Subperiod panels (2000–09 / 2010–19 / 2020–25)
# - ES95 metric; data manifest
# - SAFE turnover distribution plot (histogram/spike fallback; no SciPy KDE)
# -------------------------------------------------------------------------
import os, glob, warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------- Config ----------------
DATA_DIR = "./data2"
OUT_DIR = "./out1"
K = 3
ROLL_PC_WIN = 63 # window for PCA feature (covariance)
ROLL_FIT_N = 500 # HMM fit window (days)
REFIT_EVERY = 5
RIDGE = 1e-6
CLIP = 0.5
SEED = 7
# Reporting / fairness
TARGET_VOL = 0.10 # annualized target vol for matched-vol table
COST_BPS_GRID = [5, 10, 25, 50]
IVP_WIN = 63 # lookback for inverse-vol baseline
IVP_REBAL = 21 # rebalance every N trading days
try:
from hmmlearn.hmm import GaussianHMM
except Exception as e:
raise SystemExit("pip install hmmlearn\nOriginal error: %r" % e)
np.random.seed(SEED)
os.makedirs(OUT_DIR, exist_ok=True)
# --------------- Helpers ----------------
def load_data_from_dir(data_dir: str) -> pd.DataFrame:
files = sorted(glob.glob(os.path.join(data_dir, "*.csv")))
if not files:
raise FileNotFoundError(f"No CSVs found in {data_dir}")
frames = []
for f in files:
tkr = os.path.splitext(os.path.basename(f))[0].upper()
df = pd.read_csv(f)
if "Date" not in df.columns and "date" in df.columns:
df.rename(columns={"date": "Date"}, inplace=True)
close_col = None
for c in ["Close", "Adj Close", "close", "adj_close", "AdjClose"]:
if c in df.columns:
close_col = c; break
if close_col is None and "Zamkniecie" in df.columns:
close_col = "Zamkniecie"
if close_col is None:
print(f"Skipping {tkr}: no Close/Adj Close column.")
continue
df["Date"] = pd.to_datetime(df["Date"])
df = df[["Date", close_col]].dropna()
df.rename(columns={close_col: tkr}, inplace=True)
frames.append(df)
if not frames:
raise ValueError("No usable files with Close prices.")
out = frames[0]
for df in frames[1:]:
out = out.merge(df, on="Date", how="outer")
out.sort_values("Date", inplace=True)
out.set_index("Date", inplace=True)
return out
def to_log_returns(prices: pd.DataFrame) -> pd.DataFrame:
return np.log(prices/prices.shift(1))
def mean_var_weights(mu: np.ndarray, Sigma: np.ndarray, ridge=1e-6, clip=0.5):
n = len(mu)
S = Sigma + ridge*np.eye(n)
w = np.linalg.pinv(S).dot(mu.reshape(-1,1)).ravel()
w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones(n)/n
w = np.clip(w, -clip, clip)
w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones(n)/n
return w
def ann_kpis_from_returns(rets: pd.Series, freq=252):
rets = rets.dropna()
if rets.empty:
return dict(ret=np.nan, vol=np.nan, sharpe=np.nan, maxdd=np.nan, es95=np.nan)
curve = (1+rets).cumprod()
maxdd = (curve/curve.cummax()-1).min()
ann_ret = (1+rets).prod()**(freq/len(rets)) - 1
ann_vol = rets.std()*np.sqrt(freq)
sharpe = ann_ret/ann_vol if ann_vol>0 else np.nan
q = np.quantile(rets, 0.05)
es95 = -rets[rets<=q].mean() if (rets<=q).any() else np.nan
return dict(ret=ann_ret, vol=ann_vol, sharpe=sharpe, maxdd=maxdd, es95=es95)
def pca_pc1_var_feature(R: pd.DataFrame, win=63):
"""Return:
- pc1_share_series: log(PC1 variance share) indexed by the *window end* date
- last_loadings: PC1 eigenvector (unit length) from the *last* window (for plotting)
"""
from numpy.linalg import eigh
vals, idxs = [], []
last_loadings = None
last_cols = R.columns
for i in range(win-1, len(R)):
X = R.iloc[i-win+1:i+1].dropna(how="any")
if len(X) < win:
continue
C = X.cov().values
evals, evecs = eigh(C) # ascending order
lam1 = float(evals[-1])
total = float(evals.sum()) + 1e-12
pc1_share = lam1 / total
vals.append(np.log(pc1_share + 1e-12))
idxs.append(R.index[i])
last_loadings = evecs[:, -1] # associated eigenvector
last_cols = X.columns
s = pd.Series(vals, index=idxs, name="x_pc1_share_log")
return s, (last_cols, last_loadings)
def smooth_drawdown_from_returns(rets: pd.Series):
curve = (1+rets.fillna(0)).cumprod()
dd = curve/curve.cummax() - 1.0
return curve, dd
def inverse_vol_weights(R_window: pd.DataFrame):
vol = R_window.std().replace(0, np.nan)
w = 1.0/vol
w = w.replace([np.inf, -np.inf], np.nan).fillna(0.0)
if w.sum() > 0:
w = w / w.sum()
else:
w = pd.Series(np.ones(len(R_window.columns))/len(R_window.columns), index=R_window.columns)
return w.values
def turnover_series(W: pd.DataFrame):
"""One-way turnover per day: 0.5 * sum(|w_t - w_{t-1}|)."""
Wf = W.fillna(method="ffill").fillna(0.0)
dW = (Wf - Wf.shift(1)).abs().sum(axis=1) * 0.5
dW.iloc[0] = 0.0
return dW
def apply_costs(gross: pd.Series, turnover: pd.Series, cost_bps: float):
c = cost_bps/10000.0
return gross - c*turnover
def risk_match(rets: pd.Series, target_vol=0.10, freq=252):
vol = rets.std()*np.sqrt(freq)
if not np.isfinite(vol) or vol<=0:
return rets.copy(), 1.0
alpha = target_vol/vol
return rets*alpha, alpha
# --- ORIGINAL (non-causal) regime weights estimation -------------------
def regime_weights_from_labels(R: pd.DataFrame, labels: pd.Series, ridge=1e-6, clip=0.5):
N = R.shape[1]
out = {}
for k in sorted(labels.dropna().unique()):
idx = labels.index[labels==k]
if len(idx) < max(60, N*5):
w = np.ones(N)/N
else:
sub = R.loc[idx]
mu = sub.mean().values
Sigma = sub.cov().values
w = mean_var_weights(mu, Sigma, ridge=ridge, clip=clip)
out[int(k)] = w
return out
# ---------------- Load & prep ----------------
prices = load_data_from_dir(DATA_DIR).ffill().dropna(how="all")
tickers = [c for c in prices.columns if prices[c].notna().sum()>0]
prices = prices[tickers].dropna()
R = to_log_returns(prices).dropna()
N = R.shape[1]
print(f"Loaded {N} assets: {', '.join(R.columns)}")
# Data manifest (repro)
manifest = pd.DataFrame({
"Ticker": R.columns,
"FirstDate": [R.index.min()]*N,
"LastDate": [R.index.max()]*N,
"Obs": [R[c].dropna().shape[0] for c in R.columns]
})
manifest.to_csv(os.path.join(OUT_DIR, "manifest_pca_hmm.csv"), index=False)
# Feature: log(PC1 variance share)
x_pc1, (cols_last, loadings_last) = pca_pc1_var_feature(R, win=ROLL_PC_WIN)
x_df = x_pc1.to_frame()
R_aligned = R.loc[x_df.index].copy()
dates = x_df.index
# ---------------- Rolling HMM (labels) ----------------
labels = pd.Series(index=dates, dtype=int)
trans_mats = {}
model = None
last_fit = None
for i, dt in enumerate(dates):
refit = (last_fit is None) or ((i - last_fit) >= REFIT_EVERY)
if refit:
start_idx = max(0, i - ROLL_FIT_N)
fit_slice = x_df.iloc[start_idx:i+1].values
if len(fit_slice) < max(120, ROLL_FIT_N//4):
labels.iloc[i] = 0
continue
model = GaussianHMM(n_components=K, covariance_type='diag', n_iter=200, random_state=SEED)
model.fit(fit_slice)
last_fit = i
trans_mats[dt] = model.transmat_.copy()
start_idx = max(0, i - ROLL_FIT_N)
window = x_df.iloc[start_idx:i+1].values
path = model.predict(window)
labels.iloc[i] = path[-1]
# Remap states calm→mid→stress by feature mean
means = [(k, x_df.loc[labels==k, "x_pc1_share_log"].mean()) for k in range(K)]
order = [k for (k, _) in sorted(means, key=lambda z: z[1])]
remap = {old:new for new, old in enumerate(order)}
labels = labels.map(remap)
# ---------------- Regime weights (NON-CAUSAL, original) ---------------
regime_w = regime_weights_from_labels(R_aligned, labels, ridge=RIDGE, clip=CLIP)
# ---------------- Strategy (hard regimes) ------------------------------
W_dyn = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
prev_label = None
for dt in R_aligned.index:
if prev_label is None:
W_dyn.loc[dt] = np.ones(N)/N
else:
w = regime_w.get(int(prev_label), np.ones(N)/N)
W_dyn.loc[dt] = w
prev_label = labels.loc[dt] if pd.notna(labels.loc[dt]) else prev_label
dyn_rets = (R_aligned * W_dyn.shift(1)).sum(axis=1).fillna(0.0)
dyn_turn = turnover_series(W_dyn)
# ---------------- Baselines -------------------------------------------
# Static Mean-Var
Sigma_full = R_aligned.cov().values
mu_full = R_aligned.mean().values
w_static = mean_var_weights(mu_full, Sigma_full, ridge=RIDGE, clip=CLIP)
W_static = pd.DataFrame([w_static]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
static_rets = (R_aligned @ pd.Series(w_static, index=R_aligned.columns)).fillna(0.0)
static_turn = turnover_series(W_static)
# Equal-Weight
ew_w = np.ones(N)/N
W_ew = pd.DataFrame([ew_w]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
ew_rets = (R_aligned @ pd.Series(ew_w, index=R_aligned.columns)).fillna(0.0)
ew_turn = turnover_series(W_ew)
# Constant-Vol Risk-Parity (IVP), rolling
W_ivp = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for i, dt in enumerate(R_aligned.index):
if i < IVP_WIN:
W_ivp.iloc[i] = np.ones(N)/N
else:
if i % IVP_REBAL == 0 or i == IVP_WIN:
Rw = R_aligned.iloc[i-IVP_WIN:i]
w = inverse_vol_weights(Rw)
W_ivp.iloc[i] = w
else:
W_ivp.iloc[i] = W_ivp.iloc[i-1]
W_ivp = W_ivp.fillna(method="ffill").fillna(1.0/N)
ivp_rets = (R_aligned * W_ivp.shift(1)).sum(axis=1).fillna(0.0)
ivp_turn = turnover_series(W_ivp)
# ---------------- KPIs: native (pre-cost) -----------------------------
def kpi_row(name, rets):
k = ann_kpis_from_returns(rets)
return [name, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']]
kpi_native = pd.DataFrame(
[kpi_row('Equal-Weight', ew_rets),
kpi_row('Static Mean-Var', static_rets),
kpi_row('IVP (Const-Vol RP, rolling)', ivp_rets),
kpi_row('PCA-HMM Hard-Regime', dyn_rets)],
columns=['Strategy','Ann.Return','Ann.Vol','Sharpe','MaxDD','ES95_day']
)
kpi_native_path = os.path.join(OUT_DIR, "kpi_pca_native.csv")
kpi_native.to_csv(kpi_native_path, index=False)
print("\n=== KPI (native, pre-cost) ===")
print(kpi_native.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
print(f"Saved → {kpi_native_path}")
# ---------------- Risk-matched (target vol) headline table ------------
def matched_table(name, rets, turn):
rm_rets, alpha = risk_match(rets, TARGET_VOL)
row = ann_kpis_from_returns(rm_rets)
return name, rm_rets, turn*alpha, alpha, row
rows = []
matched_series = {}
turn_series = {}
alphas = {}
for name, rets, turn in [
('Equal-Weight', ew_rets, ew_turn),
('Static Mean-Var', static_rets, static_turn),
('IVP (Const-Vol RP, rolling)', ivp_rets, ivp_turn),
('PCA-HMM Hard-Regime', dyn_rets, dyn_turn),
]:
nm, rm_rets, rm_turn, alpha, row = matched_table(name, rets, turn)
matched_series[nm] = rm_rets
turn_series[nm] = rm_turn
alphas[nm] = alpha
rows.append([nm, row['ret'], row['vol'], row['sharpe'], row['maxdd'], row['es95'], alpha])
kpi_matched = pd.DataFrame(rows, columns=['Strategy','Ann.Return','Ann.Vol','Sharpe','MaxDD','ES95_day','Alpha(scale)'])
kpi_matched_path = os.path.join(OUT_DIR, "kpi_pca_matchedVol.csv")
kpi_matched.to_csv(kpi_matched_path, index=False)
print("\n=== KPI (risk-matched @ {:.0f}%) ===".format(TARGET_VOL*100))
print(kpi_matched.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
print(f"Saved → {kpi_matched_path}")
# ---------------- Net-of-costs grids (on risk-matched) -----------------
def cost_grid_table(cost_bps_list):
out = []
for c in cost_bps_list:
for name in matched_series:
net = apply_costs(matched_series[name], turn_series[name], c)
k = ann_kpis_from_returns(net)
out.append([c, name, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
df = pd.DataFrame(out, columns=['Cost_bps','Strategy','Ann.Return','Ann.Vol','Sharpe','MaxDD','ES95_day'])
return df
kpi_costs = cost_grid_table(COST_BPS_GRID)
kpi_costs_path = os.path.join(OUT_DIR, "kpi_pca_matchedVol_costGrid.csv")
kpi_costs.to_csv(kpi_costs_path, index=False)
print(f"\nSaved net-of-costs grid → {kpi_costs_path}")
# ---------------- Turnover stats table ---------------------------------
turn_stats = []
for name, ts in turn_series.items():
turn_stats.append([name, ts.mean(), ts.median(), ts.quantile(0.9), ts.max()])
turn_df = pd.DataFrame(turn_stats, columns=['Strategy','Turnover_mean','Turnover_median','Turnover_p90','Turnover_max'])
turn_df_path = os.path.join(OUT_DIR, "turnover_stats.csv")
turn_df.to_csv(turn_df_path, index=False)
print(f"Saved turnover stats → {turn_df_path}")
# ---------------- Subperiod panels (risk-matched, 3 eras) -------------
def subperiod_slices(idx):
eras = [
("2000-01-01","2009-12-31"),
("2010-01-01","2019-12-31"),
("2020-01-01", str(idx.max().date()))
]
out = []
for s,e in eras:
sdt = pd.Timestamp(s); edt = pd.Timestamp(e)
mask = (idx>=sdt) & (idx<=edt)
if mask.any():
out.append((s, e, mask))
return out
subs_rows = []
for (s,e,mask) in subperiod_slices(R_aligned.index):
for name in matched_series:
rets = matched_series[name].loc[mask]
k = ann_kpis_from_returns(rets)
subs_rows.append([f"{s}–{e}", name, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
sub_df = pd.DataFrame(subs_rows, columns=['Period','Strategy','Ann.Return','Ann.Vol','Sharpe','MaxDD','ES95_day'])
sub_df_path = os.path.join(OUT_DIR, "kpi_pca_matchedVol_subperiods.csv")
sub_df.to_csv(sub_df_path, index=False)
print(f"Saved subperiod KPIs → {sub_df_path}")
# ---------------- Figures -------------------
# SAFE helper: turnover distribution without SciPy KDE
def plot_turnover_distribution(turn_map, out_path):
fig, ax = plt.subplots(figsize=(10,4))
for name, ts in turn_map.items():
s = pd.Series(ts).replace([np.inf, -np.inf], np.nan).dropna()
s = s[np.isfinite(s)]
if (len(s) == 0) or (s.std() <= 1e-12) or (s.nunique() < 3):
val = float(s.iloc[0]) if len(s) else 0.0
ax.axvline(val, linestyle="--", linewidth=1.5, label=f"{name} (spike)")
else:
counts, bins = np.histogram(s, bins=40, density=True)
centers = 0.5*(bins[1:] + bins[:-1])
ax.plot(centers, counts, linewidth=1.5, label=name)
ax.set_title("Turnover distribution (one-way, risk-matched)")
ax.grid(True); ax.legend()
fig.tight_layout(); fig.savefig(out_path, dpi=160)
# 0) Turnover distribution (risk-matched scaled)
plot_turnover_distribution(turn_series, os.path.join(OUT_DIR, "fig_turnover_density.png"))
# 1) PC1 variance share (level, not log)
pc1_share_level = np.exp(x_pc1) # back from log
plt.figure(figsize=(12,3.2))
plt.plot(pc1_share_level.index, pc1_share_level.values)
plt.title("PC1 Variance Share (rolling covariance)")
plt.grid(True); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, "fig_pc1_share.png"), dpi=160)
# 2) EW price w/ regime ribbon
ew_price_curve = (1 + (R_aligned.mean(axis=1).fillna(0.0))).cumprod()
fig, ax = plt.subplots(figsize=(12,4))
ax.plot(ew_price_curve.index, ew_price_curve.values, linewidth=1.2)
ax.set_title("EW Price with Regime Ribbon (PCA-HMM, Hard)")
ax.grid(True)
reg_series = labels.reindex(ew_price_curve.index).fillna(method='ffill')
last = None
for i, dt in enumerate(ew_price_curve.index):
r = int(reg_series.iloc[i]) if not pd.isna(reg_series.iloc[i]) else None
if i==0:
last = (dt, r)
elif r != last[1]:
ax.axvspan(last[0], ew_price_curve.index[i-1], alpha=0.10)
last = (dt, r)
ax.axvspan(last[0], ew_price_curve.index[-1], alpha=0.10)
fig.tight_layout()
fig.savefig(os.path.join(OUT_DIR, "fig_regimes_ribbon.png"), dpi=160)
# 3) Transition matrix heatmap (last fit)
if len(trans_mats):
last_A = list(trans_mats.values())[-1]
plt.figure(figsize=(4,3))
plt.imshow(last_A, aspect='auto')
plt.colorbar()
plt.title("Transition Matrix (last fit)")
plt.xticks(range(K), [f"{i}" for i in range(K)])
plt.yticks(range(K), [f"{i}" for i in range(K)])
plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, "fig_transmat_heat.png"), dpi=160)
# 4) Regime MV weights (bar) — display (last window)
plt.figure(figsize=(10,3.2))
disp_reg_w = {}
Lw = labels.iloc[-ROLL_FIT_N:]
Rw = R_aligned.iloc[-ROLL_FIT_N:]
for k in range(K):
idx_k = Lw.index[Lw==k]
if len(idx_k) < max(60, N*5):
disp_reg_w[k] = np.ones(N)/N
else:
mu = Rw.loc[idx_k].mean().values
Sigma = Rw.loc[idx_k].cov().values
disp_reg_w[k] = mean_var_weights(mu, Sigma, ridge=RIDGE, clip=CLIP)
bar_data = np.array([disp_reg_w.get(k, np.ones(N)/N) for k in range(K)])
xpos = np.arange(N); width = 0.8 / K
for k in range(K):
plt.bar(xpos + k*width, bar_data[k], width=width, label=f"Regime {k}")
plt.xticks(xpos + width*(K-1)/2, R_aligned.columns, rotation=0)
plt.title("Mean–Variance Weights by Regime (display, last window)")
plt.legend(ncol=min(3,K)); plt.grid(True, axis='y'); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, "fig_regime_weights_bar.png"), dpi=160)
# 5) PC1 loadings (last window)
if loadings_last is not None:
plt.figure(figsize=(10,3.2))
vals = np.abs(loadings_last)
xpos = np.arange(len(cols_last))
plt.bar(xpos, vals)
plt.xticks(xpos, cols_last, rotation=0)
plt.title("PC1 Loadings (last rolling window)")
plt.grid(True, axis='y')
plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, "fig_pca_loadings_last.png"), dpi=160)
# 6) Log cumulative wealth (risk-matched, pre-cost)
plt.figure(figsize=(10,4))
for name, rets in matched_series.items():
curve = (1+rets).cumprod()
plt.plot(curve.index, np.log(curve), label=name)
plt.title("Log Cumulative Wealth (risk-matched, pre-cost)")
plt.legend(); plt.grid(True); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, "fig_cum_wealth_matched.png"), dpi=160)
# 7) Drawdowns (risk-matched, pre-cost)
def ddplot_from_rets(name, rets):
curve, dd = smooth_drawdown_from_returns(rets)
plt.figure(figsize=(10,3))
plt.plot(dd.index, dd.values)
plt.title(f"Drawdown (risk-matched): {name}")
plt.grid(True); plt.tight_layout()
f = f"fig_drawdown_matched_{name.replace(' ','_').replace('-','').lower()}.png"
plt.savefig(os.path.join(OUT_DIR, f), dpi=160)
for name, rets in matched_series.items():
ddplot_from_rets(name, rets)
print(f"\nAll outputs saved in: {os.path.abspath(OUT_DIR)}")
Loaded 13 assets: ACWI, BIL, GLD, GOVT, IEF, LQD, QQQ, SHV, SPY, TIP, TLT, URTH, VTI
Model is not converging. Current: 412.05708234156646 is not greater than 412.25667257937704. Delta is -0.1995902378105825 Model is not converging. Current: 402.7170148008754 is not greater than 402.72015698867034. Delta is -0.0031421877949355803 Model is not converging. Current: 548.3737935063608 is not greater than 548.37978827715. Delta is -0.0059947707892433755 Model is not converging. Current: 458.04369356984495 is not greater than 458.04467758021275. Delta is -0.0009840103678016021 Model is not converging. Current: 464.83491798200885 is not greater than 464.84100816510386. Delta is -0.0060901830950115254 Model is not converging. Current: 469.7994319010304 is not greater than 469.80392097771005. Delta is -0.004489076679647042 Model is not converging. Current: 474.4007384413255 is not greater than 474.403379971688. Delta is -0.0026415303625526576 Model is not converging. Current: 479.0086608555064 is not greater than 479.0094556877677. Delta is -0.0007948322613060554 Model is not converging. Current: 493.47385179509746 is not greater than 493.474682429957. Delta is -0.0008306348595397139 Model is not converging. Current: 498.94655364230124 is not greater than 498.94681609077185. Delta is -0.00026244847060752363 Model is not converging. Current: 554.5222344478406 is not greater than 554.5231419669689. Delta is -0.0009075191283045569 Model is not converging. Current: 648.6166496695652 is not greater than 648.6230340005968. Delta is -0.006384331031654256 Model is not converging. Current: 650.6273258593816 is not greater than 650.6339801327547. Delta is -0.006654273373101205 Model is not converging. Current: 601.9490902850579 is not greater than 601.9548589034861. Delta is -0.005768618428191985 Model is not converging. Current: 613.985342592658 is not greater than 613.9919069074178. Delta is -0.006564314759884837 Model is not converging. Current: 613.7151832723503 is not greater than 613.7166058711778. Delta is -0.0014225988275029522 Model is not converging. Current: 618.8342482801512 is not greater than 618.8343369381384. Delta is -8.865798724855267e-05 Model is not converging. Current: 576.2978718078684 is not greater than 576.3065545490861. Delta is -0.008682741217626244 Model is not converging. Current: 563.0572157864834 is not greater than 563.0594077707646. Delta is -0.0021919842811257695 Model is not converging. Current: 656.1601844149233 is not greater than 656.1717713153429. Delta is -0.011586900419615631 Model is not converging. Current: 654.4719437823305 is not greater than 654.4723628116546. Delta is -0.00041902932412085647 Model is not converging. Current: 689.4799488728324 is not greater than 689.4996453698827. Delta is -0.019696497050290418 Model is not converging. Current: 674.3351453065084 is not greater than 674.3443046440252. Delta is -0.009159337516848609 Model is not converging. Current: 674.1392438012732 is not greater than 674.1427291401491. Delta is -0.0034853388758619985 Model is not converging. Current: 708.9170892738368 is not greater than 708.9331318065316. Delta is -0.016042532694768852 Model is not converging. Current: 680.8170834978882 is not greater than 680.8191496572202. Delta is -0.002066159332002826 Model is not converging. Current: 704.4710470524119 is not greater than 704.4773219052396. Delta is -0.0062748528276870275 Model is not converging. Current: 708.5815385020234 is not greater than 708.5906692437209. Delta is -0.009130741697504163 Model is not converging. Current: 714.7325666648908 is not greater than 714.7349968885936. Delta is -0.002430223702845069 Model is not converging. Current: 616.2086402925729 is not greater than 616.2150636334702. Delta is -0.006423340897299568 Model is not converging. Current: 648.9979551877134 is not greater than 648.9989369203738. Delta is -0.0009817326604206755 Model is not converging. Current: 807.6300509384655 is not greater than 807.6384437711062. Delta is -0.00839283264065216 Model is not converging. Current: 808.9197032423106 is not greater than 808.9254432464697. Delta is -0.00574000415917908 Model is not converging. Current: 817.7482139748148 is not greater than 817.7488497110867. Delta is -0.0006357362718745208 Model is not converging. Current: 842.1303118904967 is not greater than 842.1315200774062. Delta is -0.0012081869094799913 Model is not converging. Current: 855.1886766025202 is not greater than 855.2047389166711. Delta is -0.016062314150872226 Model is not converging. Current: 853.6677937560358 is not greater than 853.68758976003. Delta is -0.019796003994201783 Model is not converging. Current: 843.0689164515139 is not greater than 843.0719885288169. Delta is -0.003072077302931575 Model is not converging. Current: 843.8445164006804 is not greater than 843.845370168209. Delta is -0.0008537675286106605 Model is not converging. Current: 852.7703008049232 is not greater than 852.7743746773743. Delta is -0.004073872451158422 Model is not converging. Current: 872.2837907737804 is not greater than 872.2896560901702. Delta is -0.005865316389758846 Model is not converging. Current: 878.6129090009238 is not greater than 878.6296856306952. Delta is -0.016776629771470652 Model is not converging. Current: 910.8490765677682 is not greater than 910.8587981314727. Delta is -0.00972156370448829 Model is not converging. Current: 920.0616191782657 is not greater than 920.063827954814. Delta is -0.0022087765482865507 Model is not converging. Current: 919.095514233573 is not greater than 919.1398364074405. Delta is -0.04432217386749926 Model is not converging. Current: 917.6305829876269 is not greater than 917.688896354362. Delta is -0.05831336673509213 Model is not converging. Current: 913.5313810322899 is not greater than 913.5793731806768. Delta is -0.04799214838692478 Model is not converging. Current: 914.6570669094224 is not greater than 914.6835833063179. Delta is -0.026516396895544858 Model is not converging. Current: 907.2671241149268 is not greater than 907.2994594092937. Delta is -0.03233529436693061 Model is not converging. Current: 910.8472404058243 is not greater than 910.8726385659303. Delta is -0.025398160105964962 Model is not converging. Current: 914.7052019674232 is not greater than 914.7750731462587. Delta is -0.0698711788354558 Model is not converging. Current: 864.9538674984052 is not greater than 864.9767809476525. Delta is -0.0229134492473122 Model is not converging. Current: 741.3025152468533 is not greater than 741.311650526179. Delta is -0.009135279325732881 Model is not converging. Current: 747.0330781434175 is not greater than 747.0684088259574. Delta is -0.03533068253989313 Model is not converging. Current: 747.5537055235167 is not greater than 747.555151791229. Delta is -0.001446267712367444 Model is not converging. Current: 767.1673054551487 is not greater than 767.181068173419. Delta is -0.013762718270299956 Model is not converging. Current: 772.3067865746947 is not greater than 772.3132942798785. Delta is -0.006507705183821599 Model is not converging. Current: 783.8850194086682 is not greater than 783.927797921195. Delta is -0.04277851252675191 Model is not converging. Current: 793.266189928237 is not greater than 793.290970168251. Delta is -0.02478024001391077 Model is not converging. Current: 802.3028458280426 is not greater than 802.3146546958699. Delta is -0.011808867827312497 Model is not converging. Current: 807.1436876705543 is not greater than 807.1591765394735. Delta is -0.0154888689191921 Model is not converging. Current: 810.5752511766023 is not greater than 810.5813260425551. Delta is -0.0060748659528826465 Model is not converging. Current: 821.7273395478358 is not greater than 821.7307527961469. Delta is -0.0034132483111761758 Model is not converging. Current: 827.9244985074381 is not greater than 827.9528469281752. Delta is -0.028348420737074775 Model is not converging. Current: 835.8103198323417 is not greater than 835.8585923483879. Delta is -0.048272516046154124 Model is not converging. Current: 847.4534753481902 is not greater than 847.5026036515641. Delta is -0.04912830337389096 Model is not converging. Current: 853.4859531129838 is not greater than 853.5392008449256. Delta is -0.05324773194183763 Model is not converging. Current: 859.2298101041796 is not greater than 859.2347896297322. Delta is -0.0049795255525850735 Model is not converging. Current: 886.5173713872819 is not greater than 886.5177155472784. Delta is -0.0003441599965299247 Model is not converging. Current: 889.0659516098925 is not greater than 889.067108567774. Delta is -0.0011569578814487613 Model is not converging. Current: 891.8801274777119 is not greater than 891.8849322486119. Delta is -0.004804770900022959 Model is not converging. Current: 888.896785396754 is not greater than 888.9021897602139. Delta is -0.005404363459888373 Model is not converging. Current: 887.9299939293911 is not greater than 887.9347266108393. Delta is -0.0047326814482175905 Model is not converging. Current: 867.5176948424016 is not greater than 867.5225249923326. Delta is -0.004830149930967309 Model is not converging. Current: 872.3943207985001 is not greater than 872.4004505907086. Delta is -0.0061297922085259415 Model is not converging. Current: 872.904772488802 is not greater than 872.910368388421. Delta is -0.005595899619038391 Model is not converging. Current: 874.5897691526415 is not greater than 874.5956475562793. Delta is -0.005878403637893825 Model is not converging. Current: 900.0019490116481 is not greater than 900.004941824777. Delta is -0.0029928131289125304 Model is not converging. Current: 910.1375998882303 is not greater than 910.1472405844145. Delta is -0.009640696184192166 Model is not converging. Current: 908.8200472438695 is not greater than 908.8253864291545. Delta is -0.005339185285038184 Model is not converging. Current: 906.0096393890362 is not greater than 906.0175465931312. Delta is -0.007907204094976805 Model is not converging. Current: 905.2579314093196 is not greater than 905.2627419622497. Delta is -0.004810552930166523 Model is not converging. Current: 945.5992188371682 is not greater than 945.600751273584. Delta is -0.0015324364158004755 Model is not converging. Current: 945.6710080866096 is not greater than 945.6985849884111. Delta is -0.027576901801467102 Model is not converging. Current: 946.8299241010748 is not greater than 946.8439878446072. Delta is -0.014063743532460649 Model is not converging. Current: 937.5310918771877 is not greater than 937.5443632747306. Delta is -0.013271397542894192 Model is not converging. Current: 931.0212196207847 is not greater than 931.0416781259825. Delta is -0.02045850519778014 Model is not converging. Current: 925.8114553463162 is not greater than 925.8289009561839. Delta is -0.01744560986776378 Model is not converging. Current: 925.5098305801708 is not greater than 925.5107867944532. Delta is -0.0009562142823824615 Model is not converging. Current: 923.0202503962339 is not greater than 923.0267896254225. Delta is -0.006539229188547324 Model is not converging. Current: 919.8072935804615 is not greater than 919.8166666734962. Delta is -0.009373093034696467 Model is not converging. Current: 926.4550415142594 is not greater than 926.4618845082546. Delta is -0.0068429939951784036 Model is not converging. Current: 935.446211722793 is not greater than 935.4660563738701. Delta is -0.019844651077050912 Model is not converging. Current: 933.8331962047333 is not greater than 933.8516910556814. Delta is -0.01849485094805914 Model is not converging. Current: 931.8850770976186 is not greater than 931.9021875736622. Delta is -0.017110476043626477 Model is not converging. Current: 928.2211334548526 is not greater than 928.2345594299329. Delta is -0.013425975080281205 Model is not converging. Current: 928.1656241651715 is not greater than 928.1772466963597. Delta is -0.011622531188208995 Model is not converging. Current: 926.2519906966522 is not greater than 926.2555130090212. Delta is -0.0035223123690002467 Model is not converging. Current: 925.8051061991202 is not greater than 925.8119030426603. Delta is -0.0067968435400871385 Model is not converging. Current: 932.6100126693772 is not greater than 932.6295007565639. Delta is -0.019488087186687153 Model is not converging. Current: 936.2844897132971 is not greater than 936.3065868385582. Delta is -0.022097125261097972 Model is not converging. Current: 942.8088491847075 is not greater than 942.8305681788586. Delta is -0.021718994151115112 Model is not converging. Current: 937.5861367518366 is not greater than 937.6030378662005. Delta is -0.01690111436380448 Model is not converging. Current: 935.2943612574704 is not greater than 935.3065056812749. Delta is -0.01214442380444325 Model is not converging. Current: 928.4146546688941 is not greater than 928.4208806932622. Delta is -0.006226024368174876 Model is not converging. Current: 926.7893408786553 is not greater than 926.796599291579. Delta is -0.0072584129237611705 Model is not converging. Current: 913.9219497992952 is not greater than 913.9296641476728. Delta is -0.007714348377589886 Model is not converging. Current: 910.4012678712803 is not greater than 910.406751128743. Delta is -0.005483257462742586 Model is not converging. Current: 904.7339778157741 is not greater than 904.7369240473801. Delta is -0.0029462316059607474 Model is not converging. Current: 895.6552862172606 is not greater than 895.6626076718456. Delta is -0.007321454585053289 Model is not converging. Current: 890.3833132029898 is not greater than 890.3920564026263. Delta is -0.008743199636455756 Model is not converging. Current: 883.6239387285055 is not greater than 883.6265558304564. Delta is -0.0026171019508183235 Model is not converging. Current: 868.4684691390399 is not greater than 868.4732295631547. Delta is -0.004760424114806483 Model is not converging. Current: 864.6202427039213 is not greater than 864.6261842668736. Delta is -0.005941562952216373 Model is not converging. Current: 867.9676120760457 is not greater than 867.9722618350156. Delta is -0.004649758969890172 Model is not converging. Current: 866.8456577379745 is not greater than 866.8502222137438. Delta is -0.004564475769257115 Model is not converging. Current: 866.5303963557111 is not greater than 866.5350120031078. Delta is -0.004615647396690292 Model is not converging. Current: 866.318337333943 is not greater than 866.3236573842102. Delta is -0.005320050267187071 Model is not converging. Current: 863.4092686713299 is not greater than 863.4154556578433. Delta is -0.006186986513398551 Model is not converging. Current: 857.035615444409 is not greater than 857.0425915121772. Delta is -0.006976067768164285 Model is not converging. Current: 849.0641831086241 is not greater than 849.0678534163704. Delta is -0.0036703077463471345 Model is not converging. Current: 846.5767774036368 is not greater than 846.6210996489799. Delta is -0.044322245343096256 Model is not converging. Current: 762.1373723644891 is not greater than 762.1469870126668. Delta is -0.009614648177716845 Model is not converging. Current: 745.1315398590408 is not greater than 745.1359356929602. Delta is -0.004395833919375036 Model is not converging. Current: 760.1161165606073 is not greater than 760.1193172246923. Delta is -0.003200664084943128
=== KPI (native, pre-cost) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day
Equal-Weight 0.0602 0.0728 0.8267 -0.2053 0.0110
Static Mean-Var 0.0182 0.0098 1.8662 -0.0132 0.0014
IVP (Const-Vol RP, rolling) 0.0199 0.0118 1.6913 -0.0216 0.0016
PCA-HMM Hard-Regime 0.0217 0.0115 1.8841 -0.0129 0.0017
Saved → ./out1/kpi_pca_native.csv
=== KPI (risk-matched @ 10%) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day Alpha(scale)
Equal-Weight 0.0821 0.1000 0.8211 -0.2728 0.0151 1.3737
Static Mean-Var 0.1977 0.1000 1.9767 -0.1307 0.0146 10.2274
IVP (Const-Vol RP, rolling) 0.1771 0.1000 1.7707 -0.1717 0.0136 8.4931
PCA-HMM Hard-Regime 0.1995 0.1000 1.9952 -0.1076 0.0144 8.6858
Saved → ./out1/kpi_pca_matchedVol.csv
Saved net-of-costs grid → ./out1/kpi_pca_matchedVol_costGrid.csv
Saved turnover stats → ./out1/turnover_stats.csv
Saved subperiod KPIs → ./out1/kpi_pca_matchedVol_subperiods.csv
All outputs saved in: /Users/jeremy.duriez/Desktop/recherche2/out1
Prob-HMM defensive core + Momentum sleeve (K=3, 500d fit, refit 3d)¶
#!/usr/bin/env python3
# Model 2 — Hybrid: Prob-HMM defensive core + Momentum sleeve (FAIR VERSION)
# -------------------------------------------------------------------------
# Your requests:
# - Keep original (non-causal) regime weight estimation (full-sample by labels) ✅
# Kept fixes:
# - Constant-vol risk-parity baseline (IVP, rolling)
# - Risk-matched (10% ann vol) headline table + net-of-costs grid (5/10/25/50 bps)
# - Subperiod KPIs (2000–09 / 2010–19 / 2020–25)
# - ES95 (daily) metric across tables; HMM calibration (Brier, ECE)
# - Data manifest for reproducibility
# - SAFE turnover distribution plot (histogram/spike fallback; no SciPy KDE) ✅
# -------------------------------------------------------------------------
import os, glob, warnings, numpy as np, pandas as pd, matplotlib.pyplot as plt
from datetime import datetime
warnings.filterwarnings("ignore")
# ---------------- Config
DATA_DIR="./data2"; OUT_DIR="./out2"
K=3; ROLL_VOL_N=21; ROLL_FIT_N=500; REFIT_EVERY=3
RIDGE=1e-6; CLIP=0.5; SEED=7
MOMO_LB=126; MOMO_MIN_OBS=90; SOFTMAX_TAU=3.0; IVOL_WIN=21
GAMMA_MAX=0.40
# Fairness / reporting
TARGET_VOL = 0.10 # risk-matched target vol (annualized)
COST_BPS_GRID = [5,10,25,50]
IVP_WIN = 63 # inverse-vol lookback
IVP_REBAL = 21 # rebalance every N days
np.random.seed(SEED); os.makedirs(OUT_DIR, exist_ok=True)
from hmmlearn.hmm import GaussianHMM
# ---------------- IO / utils
def load_data_from_dir(path):
files=sorted(glob.glob(os.path.join(path,"*.csv"))); frames=[]
for f in files:
tkr=os.path.splitext(os.path.basename(f))[0].upper()
df=pd.read_csv(f)
if "Date" not in df.columns and "date" in df.columns: df=df.rename(columns={"date":"Date"})
close_col=None
for c in ["Close","Adj Close","close","adj_close","AdjClose"]:
if c in df.columns: close_col=c; break
if close_col is None: continue
frames.append(df[["Date",close_col]].rename(columns={"Date":"Date", close_col:tkr}))
if not frames: raise FileNotFoundError("No usable CSVs.")
out=frames[0]
for df in frames[1:]: out=out.merge(df,on="Date",how="outer")
out["Date"]=pd.to_datetime(out["Date"]); out=out.sort_values("Date").set_index("Date").ffill().dropna(how="all")
return out
def to_log_returns(P): return np.log(P/P.shift(1))
def ew_vol_feature(R, n=21):
ew=R.mean(axis=1); vol=ew.rolling(n).std()
return np.log(vol + 1e-12).rename("x_ewvol")
def mean_var_weights(mu, Sigma, ridge=1e-6, clip=0.5):
S=Sigma + ridge*np.eye(len(mu))
w=np.linalg.pinv(S) @ mu.reshape(-1,1); w=w.ravel()
w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones_like(w)/len(w)
w = np.clip(w, -clip, clip); w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones_like(w)/len(w)
return w
def ann_kpis_from_returns(rets, freq=252):
rets=rets.dropna()
if rets.empty: return dict(ret=np.nan, vol=np.nan, sharpe=np.nan, maxdd=np.nan, es95=np.nan)
curve=(1+rets).cumprod()
maxdd=(curve/curve.cummax()-1).min()
ann_ret=(1+rets).prod()**(freq/len(rets)) - 1
ann_vol=rets.std()*np.sqrt(freq)
sharpe=ann_ret/ann_vol if ann_vol>0 else np.nan
q=np.quantile(rets, 0.05)
es95 = -rets[rets<=q].mean() if (rets<=q).any() else np.nan
return dict(ret=ann_ret, vol=ann_vol, sharpe=sharpe, maxdd=maxdd, es95=es95)
def smooth_drawdown_from_returns(rets):
curve=(1+rets.fillna(0)).cumprod()
dd=curve/curve.cummax()-1.0
return curve, dd
def softmax(x, tau=1.0):
x=np.asarray(x,dtype=float); x=x-np.nanmax(x); ex=np.exp(x/max(1e-12,tau))
ex[np.isnan(ex)]=0.0; s=ex.sum(); return ex/s if s>0 else np.zeros_like(ex)
def inverse_vol_weights(R_window: pd.DataFrame):
vol = R_window.std().replace(0, np.nan)
w = 1.0/vol
w = w.replace([np.inf,-np.inf], np.nan).fillna(0.0)
if w.sum()>0: w = w/w.sum()
else: w = pd.Series(np.ones(len(R_window.columns))/len(R_window.columns), index=R_window.columns)
return w.values
def turnover_series(W: pd.DataFrame):
Wf=W.fillna(method="ffill").fillna(0.0)
dW=(Wf - Wf.shift(1)).abs().sum(axis=1)*0.5
dW.iloc[0]=0.0
return dW
def apply_costs(gross: pd.Series, turnover: pd.Series, cost_bps: float):
c=cost_bps/10000.0
return gross - c*turnover
def risk_match(rets: pd.Series, target_vol=0.10, freq=252):
vol = rets.std()*np.sqrt(freq)
if not np.isfinite(vol) or vol<=0: return rets.copy(), 1.0
alpha = target_vol/vol
return rets*alpha, alpha
# ORIGINAL (non-causal) per-regime MV weights from full-sample labels
def regime_weights_from_labels(R, labels, ridge=1e-6, clip=0.5):
out={}; N=R.shape[1]
for k in sorted(labels.dropna().unique()):
idx=labels.index[labels==k]
if len(idx)<max(60,N*5): out[int(k)]=np.ones(N)/N
else:
sub=R.loc[idx]; mu=sub.mean().values; Sigma=sub.cov().values
out[int(k)] = mean_var_weights(mu,Sigma,ridge,clip)
return out
# ---------------- Data & feature
prices=load_data_from_dir(DATA_DIR)
R=to_log_returns(prices).dropna()
N=R.shape[1]; tickers=list(R.columns)
# Manifest (repro)
manifest = pd.DataFrame({
"Ticker": R.columns,
"FirstDate": [R.index.min()]*N,
"LastDate": [R.index.max()]*N,
"Obs": [R[c].dropna().shape[0] for c in R.columns]
})
os.makedirs(OUT_DIR, exist_ok=True)
manifest.to_csv(os.path.join(OUT_DIR,"m2_manifest.csv"), index=False)
x=ew_vol_feature(R, ROLL_VOL_N).dropna()
x_df=x.to_frame(); R_aligned=R.loc[x_df.index]; dates=x_df.index
print(f"Loaded {N} assets, {len(dates)} obs.")
# ---------------- Prob-HMM core (filtered & one-step-ahead)
labels=pd.Series(index=dates, dtype=int)
filt_prob=pd.DataFrame(index=dates, columns=[f"p{k}" for k in range(K)], dtype=float)
ahead_prob=filt_prob.copy()
model=None; last_fit=None
for i, dt in enumerate(dates):
refit=(last_fit is None) or ((i-last_fit)>=REFIT_EVERY)
if refit:
s=max(0,i-ROLL_FIT_N); X=x_df.iloc[s:i+1].values
if len(X)<max(120, ROLL_FIT_N//4): continue
model=GaussianHMM(n_components=K, covariance_type='diag', n_iter=200, random_state=SEED); model.fit(X); last_fit=i
s=max(0,i-ROLL_FIT_N); X=x_df.iloc[s:i+1].values
path=model.predict(X); post=model.predict_proba(X)
labels.iloc[i]=path[-1]; alpha=post[-1]; filt_prob.iloc[i]=alpha; ahead_prob.iloc[i]=alpha @ model.transmat_
# ---------------- Remap calm→mid→stress (by feature mean)
means=[(k, x_df.loc[labels==k,"x_ewvol"].mean()) for k in range(K)]
order=[k for (k,_) in sorted(means, key=lambda z:z[1])]
remap={old:new for new,old in enumerate(order)}
labels=labels.map(remap)
filt_prob=filt_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
ahead_prob=ahead_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
# ---------------- Defensive core weights (NON-CAUSAL, original)
regime_w = regime_weights_from_labels(R_aligned, labels, RIDGE, CLIP)
W_core=pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for dt in W_core.index:
p=ahead_prob.loc[dt]
if p.isna().any(): W_core.loc[dt]=np.ones(N)/N; continue
w=np.zeros(N)
for k in range(K): w += float(p[f"p{k}"])*regime_w.get(k, np.ones(N)/N)
w=np.clip(w, -CLIP, CLIP); ssum=np.abs(w).sum(); W_core.loc[dt]= w/(ssum if ssum>0 else 1.0)
# ---------------- Momentum sleeve (unchanged logic)
rets=R_aligned.copy(); price_proxy=(1+rets).cumprod()
roll_vol=rets.rolling(IVOL_WIN).std()
roll_ret=price_proxy/price_proxy.shift(MOMO_LB) - 1.0
momo_score=(roll_ret/(roll_vol+1e-12)).replace([np.inf,-np.inf], np.nan).fillna(0.0)
momo_score=momo_score.sub(momo_score.mean(axis=1), axis=0)
momo_score=momo_score.div(momo_score.std(axis=1)+1e-12, axis=0)
inv_vol=1.0/(roll_vol+1e-12); inv_vol=inv_vol.div(inv_vol.sum(axis=1), axis=0).fillna(0.0)
W_sleeve=pd.DataFrame(index=rets.index, columns=rets.columns, dtype=float)
for dt in rets.index:
srow=momo_score.loc[dt] if dt in momo_score.index else None
vrow=inv_vol.loc[dt] if dt in inv_vol.index else None
if srow is None or vrow is None or srow.isna().all() or vrow.isna().all():
W_sleeve.loc[dt]=np.ones(N)/N; continue
w_raw=softmax(srow.values, tau=SOFTMAX_TAU)
w=w_raw * vrow.values
w=w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones(N)/N
W_sleeve.loc[dt]=w
W_sleeve=W_sleeve.loc[W_core.index].fillna(method="ffill")
# ---------------- Regime-aware blend
p_stress=ahead_prob[f"p{K-1}"].reindex(W_core.index).fillna(method="ffill").clip(0,1)
gamma_t=GAMMA_MAX*(1 - p_stress)
W_final=pd.DataFrame(index=W_core.index, columns=W_core.columns, dtype=float)
for dt in W_core.index:
g=float(gamma_t.loc[dt]) if dt in gamma_t.index else 0.0
w=(1-g)*W_core.loc[dt].values + g*W_sleeve.loc[dt].values
w=np.clip(w, -CLIP, CLIP); ssum=np.abs(w).sum(); W_final.loc[dt]= w/(ssum if ssum>0 else 1.0)
# ---------------- Daily returns & turnover
ret_core =(R_aligned * W_core.shift(1)).sum(axis=1).fillna(0.0)
ret_sleeve=(R_aligned * W_sleeve.shift(1)).sum(axis=1).fillna(0.0)
ret_final =(R_aligned * W_final.shift(1)).sum(axis=1).fillna(0.0)
turn_core = turnover_series(W_core)
turn_sleeve = turnover_series(W_sleeve)
turn_final = turnover_series(W_final)
# ---------------- Baselines
# Static MV
Sigma=R_aligned.cov().values; mu=R_aligned.mean().values
w_sta=mean_var_weights(mu, Sigma, ridge=RIDGE, clip=CLIP)
W_sta = pd.DataFrame([w_sta]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
ret_sta=(R_aligned @ pd.Series(w_sta, index=R_aligned.columns)).fillna(0.0)
turn_sta = turnover_series(W_sta)
# Equal-Weight
w_ew=np.ones(N)/N
W_ew = pd.DataFrame([w_ew]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
ret_ew=(R_aligned @ pd.Series(w_ew, index=R_aligned.columns)).fillna(0.0)
turn_ew = turnover_series(W_ew)
# Constant-Vol Risk-Parity (IVP), rolling
W_ivp = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for i, dt in enumerate(R_aligned.index):
if i<IVP_WIN:
W_ivp.iloc[i]=np.ones(N)/N
else:
if i % IVP_REBAL == 0 or i == IVP_WIN:
Rw=R_aligned.iloc[i-IVP_WIN:i]
W_ivp.iloc[i]=inverse_vol_weights(Rw)
else:
W_ivp.iloc[i]=W_ivp.iloc[i-1]
W_ivp=W_ivp.fillna(method="ffill").fillna(1.0/N)
ret_ivp=(R_aligned * W_ivp.shift(1)).sum(axis=1).fillna(0.0)
turn_ivp = turnover_series(W_ivp)
# ---------------- KPIs (native, pre-cost)
def krow(name, rets):
k=ann_kpis_from_returns(rets)
return [name, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']]
kpi_native = pd.DataFrame(
[krow("Equal-Weight", ret_ew),
krow("Static MV", ret_sta),
krow("IVP (Const-Vol RP, rolling)", ret_ivp),
krow("HMM Core (prob)", ret_core),
krow("Momentum Sleeve", ret_sleeve),
krow("Hybrid Ensemble", ret_final)],
columns=["Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"]
)
kpi_native.to_csv(os.path.join(OUT_DIR,"m2_kpi_native.csv"), index=False)
print("\n=== Model 2 KPIs (native, pre-cost) ===")
print(kpi_native.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
# ---------------- Risk-matched table (10% vol)
def matched_series_and_alpha(name, rets, turn):
rm, alpha = risk_match(rets, TARGET_VOL)
return name, rm, turn*alpha, alpha, ann_kpis_from_returns(rm)
matched = {}
turn_m = {}
alphas = {}
rows=[]
for name, rets, turn in [
("Equal-Weight", ret_ew, turn_ew),
("Static MV", ret_sta, turn_sta),
("IVP (Const-Vol RP, rolling)", ret_ivp, turn_ivp),
("HMM Core (prob)", ret_core, turn_core),
("Momentum Sleeve", ret_sleeve, turn_sleeve),
("Hybrid Ensemble", ret_final, turn_final),
]:
nm, rm, tm, a, k = matched_series_and_alpha(name, rets, turn)
matched[nm]=rm; turn_m[nm]=tm; alphas[nm]=a
rows.append([nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95'], a])
kpi_matched = pd.DataFrame(rows, columns=["Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day","Alpha(scale)"])
kpi_matched.to_csv(os.path.join(OUT_DIR,"m2_kpi_matchedVol.csv"), index=False)
print("\n=== Model 2 KPIs (risk-matched @ {:.0f}%) ===".format(TARGET_VOL*100))
print(kpi_matched.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
# ---------------- Net-of-costs grid on risk-matched
def cost_grid(cost_list):
out=[]
for c in cost_list:
for nm in matched:
net = apply_costs(matched[nm], turn_m[nm], c)
k = ann_kpis_from_returns(net)
out.append([c, nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
return pd.DataFrame(out, columns=["Cost_bps","Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"])
kpi_costs = cost_grid(COST_BPS_GRID)
kpi_costs.to_csv(os.path.join(OUT_DIR,"m2_kpi_matchedVol_costGrid.csv"), index=False)
print("\nSaved net-of-costs grid → ./out/m2_kpi_matchedVol_costGrid.csv")
# ---------------- Turnover stats
turn_stats=[]
for nm, ts in turn_m.items():
turn_stats.append([nm, ts.mean(), ts.median(), ts.quantile(0.9), ts.max()])
turn_df = pd.DataFrame(turn_stats, columns=["Strategy","Turnover_mean","Turnover_median","Turnover_p90","Turnover_max"])
turn_df.to_csv(os.path.join(OUT_DIR,"m2_turnover_stats.csv"), index=False)
# ---------------- Subperiod panels (risk-matched)
def subperiod_slices(idx):
eras=[("2000-01-01","2009-12-31"), ("2010-01-01","2019-12-31"), ("2020-01-01", str(idx.max().date()))]
out=[]
for s,e in eras:
sdt=pd.Timestamp(s); edt=pd.Timestamp(e)
mask=(idx>=sdt)&(idx<=edt)
if mask.any(): out.append((s,e,mask))
return out
subs=[]
for (s,e,mask) in subperiod_slices(R_aligned.index):
for nm in matched:
k=ann_kpis_from_returns(matched[nm].loc[mask])
subs.append([f"{s}–{e}", nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
sub_df=pd.DataFrame(subs, columns=["Period","Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"])
sub_df.to_csv(os.path.join(OUT_DIR,"m2_kpi_matchedVol_subperiods.csv"), index=False)
# ---------------- HMM calibration (one-step-ahead)
# Compare ahead_prob(t) with realized label at t+1 (shift by 1)
probs = ahead_prob.copy().dropna()
y_true = labels.reindex(probs.index).shift(-1).dropna()
probs = probs.loc[y_true.index]
# Multiclass Brier
Y = np.zeros((len(y_true), K))
for i, y in enumerate(y_true.values.astype(int)):
if 0 <= y < K: Y[i, y]=1.0
P = probs[[f"p{k}" for k in range(K)]].values
brier = np.mean(np.sum((P - Y)**2, axis=1))
# ECE (top-class, 10 bins)
conf = P.max(axis=1); pred = P.argmax(axis=1); acc = (pred==y_true.values.astype(int)).astype(float)
bins = np.linspace(0,1,11); ece=0.0
for b0, b1 in zip(bins[:-1], bins[1:]):
idx = (conf>=b0) & (conf<b1)
if idx.any():
ece += np.mean(idx) * abs(acc[idx].mean() - conf[idx].mean())
cal_df = pd.DataFrame({"Metric":["Brier","ECE_top10"], "Value":[brier, ece]})
cal_df.to_csv(os.path.join(OUT_DIR,"m2_hmm_calibration.csv"), index=False)
# ---------------- SAFE turnover distribution (no SciPy KDE)
def plot_turnover_distribution(turn_map, out_path):
fig, ax = plt.subplots(figsize=(10,4))
for name, ts in turn_map.items():
s = pd.Series(ts).replace([np.inf, -np.inf], np.nan).dropna()
s = s[np.isfinite(s)]
if (len(s) == 0) or (s.std() <= 1e-12) or (s.nunique() < 3):
val = float(s.iloc[0]) if len(s) else 0.0
ax.axvline(val, linestyle="--", linewidth=1.5, label=f"{name} (spike)")
else:
counts, bins = np.histogram(s, bins=40, density=True)
centers = 0.5*(bins[1:] + bins[:-1])
ax.plot(centers, counts, linewidth=1.5, label=name)
ax.set_title("Turnover distribution (one-way, risk-matched)")
ax.grid(True); ax.legend()
fig.tight_layout(); fig.savefig(out_path, dpi=180)
plot_turnover_distribution(turn_m, os.path.join(OUT_DIR,"m2_fig_turnover_density.png"))
# ---------------- Figures (standardized, comparable)
# Risk-matched log cumulative wealth (pre-cost)
plt.figure(figsize=(10,4))
for nm, sr in matched.items():
curve=(1+sr).cumprod()
plt.plot(curve.index, np.log(curve), label=nm)
plt.title("Log Cumulative Wealth — Risk-matched (10% vol)"); plt.legend(); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m2_fig_cum_wealth_matched.png"), dpi=180)
# Risk-matched drawdowns
def ddplot_from_rets(name, rets):
curve, dd = smooth_drawdown_from_returns(rets)
plt.figure(figsize=(10,3)); plt.plot(dd)
plt.title(f"Drawdown (risk-matched) — {name}"); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, f"m2_fig_drawdown_matched_{name.replace(' ','_').replace('-','')}.png"), dpi=180)
for nm, sr in matched.items():
ddplot_from_rets(nm, sr)
# Native wealth (pre-cost)
curve_core =(1+ret_core).cumprod()
curve_sleeve=(1+ret_sleeve).cumprod()
curve_final =(1+ret_final).cumprod()
curve_sta =(1+ret_sta).cumprod()
curve_ew =(1+ret_ew).cumprod()
plt.figure(figsize=(10,4))
plt.plot(np.log(curve_final), label="Hybrid Ensemble")
plt.plot(np.log(curve_core), label="HMM Core (prob)")
plt.plot(np.log(curve_sleeve), label="Momentum Sleeve")
plt.plot(np.log(curve_sta), label="Static MV")
plt.plot(np.log(curve_ew), label="Equal-Weight")
plt.title("Log Cumulative Wealth — Model 2 (native)"); plt.legend(); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m2_cum_wealth_native.png"), dpi=180)
# Prob ribbon
cols=[f"p{k}" for k in range(K)]
pdf=ahead_prob[cols].dropna(how="any")
plt.figure(figsize=(12,3.6))
plt.stackplot(pdf.index, pdf.values.T, labels=cols, alpha=0.9)
plt.ylim(0,1); plt.title("One-step-ahead Probabilities — Model 2"); plt.legend(ncol=K, fontsize=8, frameon=False)
plt.grid(True, alpha=.2); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m2_prob_ribbon.png"), dpi=180)
# Weights timelines
plt.figure(figsize=(11,4))
plt.stackplot(W_core.index, W_core.fillna(method="ffill").values.T, labels=tickers)
plt.title("Portfolio Weights — HMM Core (prob)"); plt.legend(loc="upper left", ncol=min(N,4), fontsize=8)
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR,"m2_weights_core.png"), dpi=180)
plt.figure(figsize=(11,4))
plt.stackplot(W_sleeve.index, W_sleeve.fillna(method="ffill").values.T, labels=tickers)
plt.title("Portfolio Weights — Momentum Sleeve"); plt.legend(loc="upper left", ncol=min(N,4), fontsize=8)
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR,"m2_weights_sleeve.png"), dpi=180)
plt.figure(figsize=(11,4))
plt.stackplot(W_final.index, W_final.fillna(method="ffill").values.T, labels=tickers)
plt.title("Portfolio Weights — Hybrid Ensemble"); plt.legend(loc="upper left", ncol=min(N,4), fontsize=8)
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR,"m2_weights_final.png"), dpi=180)
# Sleeve intensity
plt.figure(figsize=(10,3))
plt.plot(gamma_t.index, gamma_t.values)
plt.title(r"Sleeve Intensity $\gamma_t$ = GAMMA_MAX · (1 − P_{stress})"); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m2_sleeve_intensity.png"), dpi=180)
# Momentum heatmap (z-scores)
momo_sub=momo_score.loc[W_core.index].copy(); Z=momo_sub.values.T
fig, ax = plt.subplots(figsize=(12, 3 + 0.12*N))
im=ax.imshow(Z, aspect="auto", cmap="viridis", vmin=-2.5, vmax=2.5)
ax.set_yticks(range(N)); ax.set_yticklabels(tickers)
xt=np.linspace(0, Z.shape[1]-1, 8, dtype=int)
ax.set_xticks(xt); ax.set_xticklabels([momo_sub.index[i].strftime("%Y") for i in xt])
fig.colorbar(im, ax=ax, fraction=0.025, pad=0.02).set_label("Momentum z-score")
ax.set_title("Momentum Heatmap — Model 2"); fig.tight_layout()
fig.savefig(os.path.join(OUT_DIR,"m2_momo_heatmap.png"), dpi=180)
print(f"\nFigures & tables saved in: {os.path.abspath(OUT_DIR)}")
Loaded 13 assets, 3404 obs.
Model is not converging. Current: 354.1968677568538 is not greater than 354.1970358707842. Delta is -0.00016811393038551614 Model is not converging. Current: 354.429513131319 is not greater than 354.43020147967803. Delta is -0.0006883483590058859 Model is not converging. Current: 235.14949633390336 is not greater than 235.14983045425902. Delta is -0.0003341203556601613 Model is not converging. Current: 241.04640901288684 is not greater than 241.0471493121839. Delta is -0.0007402992970639843 Model is not converging. Current: 219.34526619659013 is not greater than 219.34556335749306. Delta is -0.0002971609029316369 Model is not converging. Current: 214.49601436229398 is not greater than 214.49621822328905. Delta is -0.00020386099507163635 Model is not converging. Current: 249.09560366422386 is not greater than 249.09585198536337. Delta is -0.0002483211395087892 Model is not converging. Current: 233.03566808977445 is not greater than 233.03633397871312. Delta is -0.0006658889386699229 Model is not converging. Current: 12.730394609057656 is not greater than 12.7305804804049. Delta is -0.0001858713472433493 Model is not converging. Current: -4.112930001313416 is not greater than -4.112706239098814. Delta is -0.0002237622146017415 Model is not converging. Current: -4.287934729930699 is not greater than -4.287668253373226. Delta is -0.00026647655747247256 Model is not converging. Current: -5.089796482525532 is not greater than -5.089565140890337. Delta is -0.00023134163519511475 Model is not converging. Current: -6.456079962763253 is not greater than -6.455907125397283. Delta is -0.0001728373659704019 Model is not converging. Current: -7.8114834109455025 is not greater than -7.811381662566985. Delta is -0.00010174837851728569 Model is not converging. Current: 17.548726420763575 is not greater than 17.549236372431018. Delta is -0.0005099516674427207 Model is not converging. Current: 78.56205558779848 is not greater than 78.56212918975174. Delta is -7.360195326100438e-05 Model is not converging. Current: 166.48253361906458 is not greater than 166.48298203453913. Delta is -0.00044841547455121145 Model is not converging. Current: 172.88710270855572 is not greater than 172.8881344222597. Delta is -0.0010317137039805857 Model is not converging. Current: 160.4287719262974 is not greater than 160.43127858962117. Delta is -0.002506663323771363
=== Model 2 KPIs (native, pre-cost) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day
Equal-Weight 0.0582 0.0725 0.8030 -0.2053 0.0109
Static MV 0.0179 0.0096 1.8577 -0.0135 0.0014
IVP (Const-Vol RP, rolling) 0.0171 0.0120 1.4249 -0.0306 0.0017
HMM Core (prob) 0.0239 0.0135 1.7693 -0.0306 0.0018
Momentum Sleeve 0.0162 0.0087 1.8514 -0.0184 0.0012
Hybrid Ensemble 0.0198 0.0087 2.2634 -0.0119 0.0012
=== Model 2 KPIs (risk-matched @ 10%) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day Alpha(scale)
Equal-Weight 0.0797 0.1000 0.7968 -0.2738 0.0150 1.3790
Static MV 0.1967 0.1000 1.9669 -0.1354 0.0144 10.3920
IVP (Const-Vol RP, rolling) 0.1467 0.1000 1.4666 -0.2306 0.0142 8.3263
HMM Core (prob) 0.1859 0.1000 1.8590 -0.2077 0.0135 7.4043
Momentum Sleeve 0.1961 0.1000 1.9607 -0.1948 0.0140 11.4303
Hybrid Ensemble 0.2454 0.1000 2.4544 -0.1332 0.0141 11.4413
Saved net-of-costs grid → ./out/m2_kpi_matchedVol_costGrid.csv
Figures & tables saved in: /Users/jeremy.duriez/Desktop/recherche2/out2
Model 3 — RA-HRP: Prob-HMM on features + per-regime HRP allocation¶
#!/usr/bin/env python3
# Model 3 — RA-HRP: Prob-HMM on features + per-regime HRP allocation (FAIR VERSION)
# -------------------------------------------------------------------------
# Kept core:
# - Per-regime HRP weights blended by one-step-ahead HMM probs (your original)
# Additions:
# - HRP on mixture covariance Σ̄_t = Σ_s p_s(t) (exact alternative to prob-blend)
# - Constant-vol Risk-Parity (IVP) baseline (rolling)
# - Risk-matched (10% ann vol) headline table + net-of-costs grid (5/10/25/50 bps)
# - Subperiod KPIs (2000–09 / 2010–19 / 2020–*)
# - ES95 (daily) metric, turnover stats, safe turnover plot (no SciPy KDE)
# - HMM calibration (Brier, ECE) and data manifest
# FIXES:
# - "Static HRP" now FAIR (expanding window + monthly rebalance + shrinkage; no look-ahead)
# - Portfolio returns computed as simple returns from weighted log-returns
# -------------------------------------------------------------------------
import os, glob, warnings, numpy as np, pandas as pd, matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
# ---------- Config
DATA_DIR="./data2"; OUT_DIR="./out3"
K=3; ROLL_VOL_N=21; ROLL_FIT_N=500; REFIT_EVERY=3
CLIP=0.5; SEED=7
TARGET_VOL = 0.10 # risk-matched target annual vol
COST_BPS_GRID = [5,10,25,50] # bps per one-way turnover
IVP_WIN=63; IVP_REBAL=21 # IVP baseline params
LABEL_HRP = "Static HRP" # keep legacy name for comparability
np.random.seed(SEED); os.makedirs(OUT_DIR, exist_ok=True)
from hmmlearn.hmm import GaussianHMM
# --- Shrinkage (optional, with fallback) ---
try:
from sklearn.covariance import LedoitWolf
except Exception:
LedoitWolf = None
# ---------- IO / core utils
def load_data_from_dir(path):
files=sorted(glob.glob(os.path.join(path,"*.csv"))); frames=[]
for f in files:
tkr=os.path.splitext(os.path.basename(f))[0].upper()
df=pd.read_csv(f)
if "Date" not in df.columns and "date" in df.columns: df=df.rename(columns={"date":"Date"})
for c in ["Close","Adj Close","close","adj_close","AdjClose"]:
if c in df.columns:
frames.append(df[["Date",c]].rename(columns={"Date":"Date", c:tkr})); break
if not frames: raise FileNotFoundError("No usable CSVs.")
out=frames[0]
for df in frames[1:]: out=out.merge(df,on="Date",how="outer")
out["Date"]=pd.to_datetime(out["Date"]); out=out.sort_values("Date").set_index("Date").ffill().dropna(how="all")
return out
def to_log_returns(P): return np.log(P/P.shift(1))
def ew_vol_feature(R, n=21):
ew=R.mean(axis=1); vol=ew.rolling(n).std()
return np.log(vol + 1e-12).rename("x_ewvol")
def pca_pc1_var_feature(R, win=63):
from numpy.linalg import eigh
vals, idxs = [], []
for i in range(win-1, len(R)):
X = R.iloc[i-win+1:i+1].dropna(how="any")
if len(X)<win: continue
C = X.cov().values
evals,_ = eigh(C); lam1=float(evals[-1]); total=float(evals.sum())+1e-12
vals.append(np.log(lam1/total + 1e-12)); idxs.append(R.index[i])
return pd.Series(vals, index=idxs, name="x_pc1")
def ann_kpis_from_returns(rets, freq=252):
rets = pd.Series(rets).dropna()
if rets.empty: return dict(ret=np.nan, vol=np.nan, sharpe=np.nan, maxdd=np.nan, es95=np.nan)
curve=(1+rets).cumprod()
maxdd=(curve/curve.cummax()-1).min()
ann_ret=(1+rets).prod()**(freq/len(rets)) - 1
ann_vol=rets.std()*np.sqrt(freq)
sharpe=ann_ret/ann_vol if ann_vol>0 else np.nan
q=np.quantile(rets, 0.05)
es95 = -rets[rets<=q].mean() if (rets<=q).any() else np.nan
return dict(ret=ann_ret, vol=ann_vol, sharpe=sharpe, maxdd=maxdd, es95=es95)
def smooth_drawdown_from_returns(rets):
curve=(1+pd.Series(rets).fillna(0)).cumprod()
dd=curve/curve.cummax()-1.0
return curve, dd
# ---------- HRP helpers
def inverse_variance_portfolio(Sigma):
d=np.diag(Sigma).clip(1e-12,None); iv=1.0/d; w=iv/iv.sum(); return w
def _corr_from_cov(Sigma):
std=np.sqrt(np.diag(Sigma)).clip(1e-12,None)
return Sigma/np.outer(std,std)
def _seriation_greedy(Corr):
N=Corr.shape[0]; remaining=list(range(N))
start=int(np.argmax(Corr.mean(axis=1))); order=[start]; remaining.remove(start)
while remaining:
last=order[-1]; nxt=max(remaining, key=lambda j: Corr[last,j])
order.append(nxt); remaining.remove(nxt)
return order
def hrp_weights(Sigma):
Corr=_corr_from_cov(Sigma); order=_seriation_greedy(Corr)
N=Sigma.shape[0]; w=np.ones(N); clusters=[order]
while clusters:
new=[]
for cl in clusters:
if len(cl)<=1: continue
m=len(cl)//2; c1, c2 = cl[:m], cl[m:]
S1=Sigma[np.ix_(c1,c1)]; S2=Sigma[np.ix_(c2,c2)]
w1=inverse_variance_portfolio(S1); v1=float(w1@S1@S1@w1) if False else float(w1@S1@w1)
w2=inverse_variance_portfolio(S2); v2=float(w2@S2@S2@w2) if False else float(w2@S2@w2)
a = 1.0 - v1/(v1+v2) if (v1+v2)>0 else 0.5
w[c1]*=a; w[c2]*=(1.0-a)
if len(c1)>1: new.append(c1)
if len(c2)>1: new.append(c2)
clusters=new
w=np.clip(w,0,None); s=np.abs(w).sum(); return w/(s if s>0 else 1.0)
def inverse_vol_weights(R_window: pd.DataFrame):
vol = R_window.std().replace(0, np.nan)
w = 1.0/vol
w = w.replace([np.inf,-np.inf], np.nan).fillna(0.0)
if w.sum()>0: w = w/w.sum()
else: w = pd.Series(np.ones(len(R_window.columns))/len(R_window.columns), index=R_window.columns)
return w.values
def turnover_series(W: pd.DataFrame):
Wf=W.fillna(method="ffill").fillna(0.0)
dW=(Wf - Wf.shift(1)).abs().sum(axis=1)*0.5
dW.iloc[0]=0.0
return dW
def apply_costs(gross: pd.Series, turnover: pd.Series, cost_bps: float):
c=cost_bps/10000.0
return gross - c*turnover
def risk_match(rets: pd.Series, target_vol=0.10, freq=252):
vol = rets.std()*np.sqrt(freq)
if not np.isfinite(vol) or vol<=0: return rets.copy(), 1.0
alpha = target_vol/vol
return rets*alpha, alpha
def shrink_cov_matrix(R_window: pd.DataFrame):
"""
Robust covariance:
- Ledoit–Wolf shrinkage when available
- PSD projection (eigenvalue floor) as a safety
"""
X = R_window.values
if LedoitWolf is not None and X.shape[0] > X.shape[1] + 2:
try:
S = LedoitWolf().fit(X).covariance_
except Exception:
S = R_window.cov().values
else:
S = R_window.cov().values
S = np.asarray(S, float)
S[~np.isfinite(S)] = 0.0
vals, vecs = np.linalg.eigh(S)
vals = np.clip(vals, 1e-8, None)
S_psd = (vecs * vals) @ vecs.T
return S_psd
# ---------- Data & features (two-feature HMM like paper)
prices=load_data_from_dir(DATA_DIR)
R=to_log_returns(prices).dropna()
N=R.shape[1]; tickers=list(R.columns)
# Manifest (repro)
manifest = pd.DataFrame({
"Ticker": R.columns,
"FirstDate": [R.index.min()]*N,
"LastDate": [R.index.max()]*N,
"Obs": [R[c].dropna().shape[0] for c in R.columns]
})
manifest.to_csv(os.path.join(OUT_DIR, "m3_manifest.csv"), index=False)
x1=ew_vol_feature(R, n=ROLL_VOL_N)
x2=pca_pc1_var_feature(R, win=ROLL_VOL_N)
x_df=pd.concat([x1,x2], axis=1).dropna()
R_aligned=R.loc[x_df.index]; dates=x_df.index
print(f"Loaded {N} assets, {len(dates)} obs, features={list(x_df.columns)}")
# ---------- Prob-HMM on features
labels=pd.Series(index=dates, dtype=int)
filt_prob=pd.DataFrame(index=dates, columns=[f"p{k}" for k in range(K)], dtype=float)
ahead_prob=filt_prob.copy()
model=None; last_fit=None
for i, dt in enumerate(dates):
refit=(last_fit is None) or ((i-last_fit)>=REFIT_EVERY)
if refit:
s=max(0,i-ROLL_FIT_N); X=x_df.iloc[s:i+1].values
if len(X)<max(120, ROLL_FIT_N//4): continue
model=GaussianHMM(n_components=K, covariance_type='diag', n_iter=200, random_state=SEED); model.fit(X); last_fit=i
s=max(0,i-ROLL_FIT_N); X=x_df.iloc[s:i+1].values
path=model.predict(X); post=model.predict_proba(X)
labels.iloc[i]=path[-1]; alpha=post[-1]; filt_prob.iloc[i]=alpha; ahead_prob.iloc[i]=alpha @ model.transmat_
# ---------- Remap calm→mid→stress
z1=(x_df["x_ewvol"]-x_df["x_ewvol"].mean())/(x_df["x_ewvol"].std()+1e-12)
z2=(x_df["x_pc1"] -x_df["x_pc1"].mean()) /(x_df["x_pc1"].std() +1e-12)
stress=z1+z2
means=[(k, stress[labels==k].mean()) for k in range(K)]
order=[k for (k,_) in sorted(means, key=lambda z:z[1])]
remap={old:new for new,old in enumerate(order)}
labels=labels.map(remap)
filt_prob=filt_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
ahead_prob=ahead_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
# ---------- Per-regime HRP weights (NON-CAUSAL, original)
def regime_hrp_weights(R, labels):
out={}
for k in sorted(labels.dropna().unique()):
idx=labels.index[labels==k]
if len(idx) < max(60, R.shape[1]*5):
out[int(k)] = np.ones(R.shape[1])/R.shape[1]; continue
Sigma = R.loc[idx].cov().values
try: w = hrp_weights(Sigma)
except Exception: w = inverse_variance_portfolio(Sigma)
w=np.clip(w, -CLIP, CLIP); s=np.abs(w).sum(); out[int(k)] = w/(s if s>0 else 1.0)
return out
regime_w=regime_hrp_weights(R_aligned, labels)
# ---------- Also precompute per-regime covariance matrices (for Σ̄)
regime_Sigma={}
for k in sorted(labels.dropna().unique()):
idx=labels.index[labels==k]
if len(idx) < max(60, R.shape[1]*5):
regime_Sigma[int(k)] = np.diag(R_aligned.var().values)
else:
regime_Sigma[int(k)] = R_aligned.loc[idx].cov().values
# ---------- Returns helper: weighted log-returns → simple returns
def port_simple_returns_from_log(R_log: pd.DataFrame, W: pd.DataFrame):
r_log = (R_log * W.shift(1)).sum(axis=1).fillna(0.0)
return np.expm1(r_log)
# ---------- Portfolios
# (A) Original RA-HRP: prob-blend of per-regime HRP weights
W_prob=pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for dt in W_prob.index:
p=ahead_prob.loc[dt]
if p.isna().any(): W_prob.loc[dt]=np.ones(N)/N; continue
w=np.zeros(N)
for k in range(K): w += float(p[f"p{k}"])*regime_w.get(k, np.ones(N)/N)
w=np.clip(w, -CLIP, CLIP); s=np.abs(w).sum(); W_prob.loc[dt]= w/(s if s>0 else 1.0)
# (B) Mixture-covariance HRP: HRP on Σ̄_t = Σ_s p_s(t)
W_mix=pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for dt in W_mix.index:
p=ahead_prob.loc[dt]
if p.isna().any(): W_mix.loc[dt]=np.ones(N)/N; continue
Sbar=np.zeros((N,N))
for k in range(K):
Sbar += float(p[f"p{k}"]) * regime_Sigma.get(k, np.diag(R_aligned.var().values))
try: w = hrp_weights(Sbar)
except Exception: w = inverse_variance_portfolio(Sbar)
w=np.clip(w, -CLIP, CLIP); s=np.abs(w).sum(); W_mix.loc[dt]= w/(s if s>0 else 1.0)
# Daily simple returns & turnover
ret_prob = port_simple_returns_from_log(R_aligned, W_prob)
ret_mix = port_simple_returns_from_log(R_aligned, W_mix)
turn_prob = turnover_series(W_prob)
turn_mix = turnover_series(W_mix)
# ---------- Baselines
# Equal-Weight (constant)
w_ew=np.ones(N)/N
W_ew=pd.DataFrame([w_ew]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
ret_ew=port_simple_returns_from_log(R_aligned, W_ew)
turn_ew=turnover_series(W_ew)
# "Static HRP" (FAIR): expanding window, monthly rebalance, shrinkage (no look-ahead)
HRP_MIN_WIN = max(252, ROLL_FIT_N) # at least ~1y or your model fit length
HRP_REBAL = 21 # ~monthly
W_hrp = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
last_w = np.ones(N)/N
for i, dt in enumerate(R_aligned.index):
if i < HRP_MIN_WIN:
W_hrp.iloc[i] = last_w
continue
if (i == HRP_MIN_WIN) or (i % HRP_REBAL == 0):
Rw = R_aligned.iloc[:i] # EXPANDING window up to t-1
Sigma_hat = shrink_cov_matrix(Rw)
try:
w = hrp_weights(Sigma_hat)
except Exception:
w = inverse_variance_portfolio(Sigma_hat)
w = np.clip(w, -CLIP, CLIP)
s = np.abs(w).sum()
last_w = w / (s if s > 0 else 1.0)
W_hrp.iloc[i] = last_w
W_hrp = W_hrp.fillna(method="ffill").fillna(1.0/N)
ret_hrp = port_simple_returns_from_log(R_aligned, W_hrp)
turn_hrp = turnover_series(W_hrp)
# Constant-Vol Risk-Parity (IVP), rolling
W_ivp = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for i, dt in enumerate(R_aligned.index):
if i<IVP_WIN:
W_ivp.iloc[i]=np.ones(N)/N
else:
if i % IVP_REBAL == 0 or i == IVP_WIN:
Rw=R_aligned.iloc[i-IVP_WIN:i]
W_ivp.iloc[i]=inverse_vol_weights(Rw)
else:
W_ivp.iloc[i]=W_ivp.iloc[i-1]
W_ivp=W_ivp.fillna(method="ffill").fillna(1.0/N)
ret_ivp=port_simple_returns_from_log(R_aligned, W_ivp)
turn_ivp=turnover_series(W_ivp)
# ---------- KPIs (native, pre-cost)
def krow(name, rets):
k=ann_kpis_from_returns(rets)
return [name, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']]
kpi_native = pd.DataFrame(
[krow("Equal-Weight", ret_ew),
krow(LABEL_HRP, ret_hrp), # keep label "Static HRP"
krow("IVP (Const-Vol RP, rolling)", ret_ivp),
krow("RA-HRP (prob-blend)", ret_prob),
krow("RA-HRP (mixture-cov HRP)", ret_mix)],
columns=["Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"]
)
kpi_native.to_csv(os.path.join(OUT_DIR,"m3_kpi_native.csv"), index=False)
print("\n=== Model 3 KPIs (native, pre-cost) ===")
print(kpi_native.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
# ---------- Risk-matched table (10% vol)
def matched_series_and_alpha(name, rets, turn):
rm, alpha = risk_match(rets, TARGET_VOL)
return name, rm, turn*alpha, alpha, ann_kpis_from_returns(rm)
matched = {}
turn_m = {}
alphas = {}
rows=[]
for name, rets, turn in [
("Equal-Weight", ret_ew, turn_ew),
(LABEL_HRP, ret_hrp, turn_hrp), # keep label "Static HRP"
("IVP (Const-Vol RP, rolling)", ret_ivp, turn_ivp),
("RA-HRP (prob-blend)", ret_prob, turn_prob),
("RA-HRP (mixture-cov HRP)", ret_mix, turn_mix),
]:
nm, rm, tm, a, k = matched_series_and_alpha(name, rets, turn)
matched[nm]=rm; turn_m[nm]=tm; alphas[nm]=a
rows.append([nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95'], a])
kpi_matched = pd.DataFrame(rows, columns=["Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day","Alpha(scale)"])
kpi_matched.to_csv(os.path.join(OUT_DIR,"m3_kpi_matchedVol.csv"), index=False)
print("\n=== Model 3 KPIs (risk-matched @ {:.0f}%) ===".format(TARGET_VOL*100))
print(kpi_matched.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
# ---------- Net-of-costs grid on risk-matched
def cost_grid(cost_list):
out=[]
for c in cost_list:
for nm in matched:
net = apply_costs(matched[nm], turn_m[nm], c)
k = ann_kpis_from_returns(net)
out.append([c, nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
return pd.DataFrame(out, columns=["Cost_bps","Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"])
kpi_costs = cost_grid(COST_BPS_GRID)
kpi_costs.to_csv(os.path.join(OUT_DIR,"m3_kpi_matchedVol_costGrid.csv"), index=False)
print("\nSaved net-of-costs grid → ./out/m3_kpi_matchedVol_costGrid.csv")
# ---------- Turnover stats
turn_stats=[]
for nm, ts in turn_m.items():
turn_stats.append([nm, ts.mean(), ts.median(), ts.quantile(0.9), ts.max()])
turn_df = pd.DataFrame(turn_stats, columns=["Strategy","Turnover_mean","Turnover_median","Turnover_p90","Turnover_max"])
turn_df.to_csv(os.path.join(OUT_DIR,"m3_turnover_stats.csv"), index=False)
# ---------- Subperiod panels (risk-matched)
def subperiod_slices(idx):
eras=[("2000-01-01","2009-12-31"), ("2010-01-01","2019-12-31"), ("2020-01-01", str(idx.max().date()))]
out=[]
for s,e in eras:
sdt=pd.Timestamp(s); edt=pd.Timestamp(e)
mask=(idx>=sdt)&(idx<=edt)
if mask.any(): out.append((s,e,mask))
return out
subs=[]
for (s,e,mask) in subperiod_slices(R_aligned.index):
for nm in matched:
k=ann_kpis_from_returns(matched[nm].loc[mask])
subs.append([f"{s}–{e}", nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
sub_df=pd.DataFrame(subs, columns=["Period","Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"])
sub_df.to_csv(os.path.join(OUT_DIR,"m3_kpi_matchedVol_subperiods.csv"), index=False)
# ---------- HMM calibration (one-step-ahead)
probs = ahead_prob.copy().dropna()
y_true = labels.reindex(probs.index).shift(-1).dropna()
probs = probs.loc[y_true.index]
Y = np.zeros((len(y_true), K))
for i, y in enumerate(y_true.values.astype(int)):
if 0 <= y < K: Y[i, y]=1.0
P = probs[[f"p{k}" for k in range(K)]].values
brier = np.mean(np.sum((P - Y)**2, axis=1))
conf = P.max(axis=1); pred = P.argmax(axis=1); acc = (pred==y_true.values.astype(int)).astype(float)
bins = np.linspace(0,1,11); ece=0.0
for b0, b1 in zip(bins[:-1], bins[1:]):
idx = (conf>=b0) & (conf<b1)
if idx.any():
ece += np.mean(idx) * abs(acc[idx].mean() - conf[idx].mean())
pd.DataFrame({"Metric":["Brier","ECE_top10"], "Value":[brier, ece]}).to_csv(os.path.join(OUT_DIR,"m3_hmm_calibration.csv"), index=False)
# ---------- SAFE turnover distribution (no SciPy KDE)
def plot_turnover_distribution(turn_map, out_path):
fig, ax = plt.subplots(figsize=(10,4))
for name, ts in turn_map.items():
s = pd.Series(ts).replace([np.inf, -np.inf], np.nan).dropna()
s = s[np.isfinite(s)]
if (len(s) == 0) or (s.std() <= 1e-12) or (s.nunique() < 3):
val = float(s.iloc[0]) if len(s) else 0.0
ax.axvline(val, linestyle="--", linewidth=1.5, label=f"{name} (spike)")
else:
counts, bins = np.histogram(s, bins=40, density=True)
centers = 0.5*(bins[1:] + bins[:-1])
ax.plot(centers, counts, linewidth=1.5, label=name)
ax.set_title("Turnover distribution (one-way, risk-matched)")
ax.grid(True); ax.legend()
fig.tight_layout(); fig.savefig(out_path, dpi=180)
plot_turnover_distribution(turn_m, os.path.join(OUT_DIR,"m3_fig_turnover_density.png"))
# ---------- Figures (standardized)
# Risk-matched log cumulative wealth (pre-cost)
plt.figure(figsize=(10,4))
for nm, sr in matched.items():
curve=(1+sr).cumprod()
plt.plot(curve.index, np.log(curve), label=nm)
plt.title("Log Cumulative Wealth — Risk-matched (10% vol)"); plt.legend(); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m3_fig_cum_wealth_matched.png"), dpi=180)
# Risk-matched drawdowns
def ddplot_from_rets(name, rets):
curve, dd = smooth_drawdown_from_returns(rets)
plt.figure(figsize=(10,3)); plt.plot(dd)
plt.title(f"Drawdown (risk-matched) — {name}"); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, f"m3_fig_drawdown_matched_{name.replace(' ','_').replace('-','')}.png"), dpi=180)
for nm, sr in matched.items():
ddplot_from_rets(nm, sr)
# Native wealth (pre-cost)
curve_prob=(1+ret_prob).cumprod(); curve_mix=(1+ret_mix).cumprod()
curve_hrp =(1+ret_hrp ).cumprod(); curve_ew=(1+ret_ew).cumprod()
plt.figure(figsize=(10,4))
plt.plot(np.log(curve_prob), label="RA-HRP (prob-blend)")
plt.plot(np.log(curve_mix), label="RA-HRP (mixture-cov HRP)")
plt.plot(np.log(curve_hrp), label=LABEL_HRP) # keep label "Static HRP"
plt.plot(np.log(curve_ew), label="Equal-Weight")
plt.title("Log Cumulative Wealth — Model 3 (native)"); plt.legend(); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m3_cum_wealth_native.png"), dpi=180)
# One-step-ahead probability ribbon
cols=[f"p{k}" for k in range(K)]
pdf=ahead_prob[cols].dropna(how="any")
plt.figure(figsize=(12,3.6))
plt.stackplot(pdf.index, pdf.values.T, labels=cols, alpha=0.9)
plt.ylim(0,1); plt.title("One-step-ahead Probabilities — Model 3"); plt.legend(ncol=K, fontsize=8, frameon=False)
plt.grid(True, alpha=.2); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"m3_prob_ribbon.png"), dpi=180)
# Weights timeline (prob-blend portfolio)
plt.figure(figsize=(11,4))
plt.stackplot(W_prob.index, W_prob.fillna(method="ffill").values.T, labels=tickers)
plt.title("Portfolio Weights — RA-HRP (prob-blend)"); plt.legend(loc="upper left", ncol=min(N,4), fontsize=8)
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR,"m3_weights_prob_blend.png"), dpi=180)
# ---------- Methods note for reproducibility (kept label "Static HRP")
pd.DataFrame({
"Strategy": [LABEL_HRP],
"Method": ["Expanding window, monthly rebalance, Ledoit–Wolf shrinkage, PSD; no look-ahead"]
}).to_csv(os.path.join(OUT_DIR, "m3_method_notes.csv"), index=False)
print(f"\nFigures & tables saved in: {os.path.abspath(OUT_DIR)}")
Model is not converging. Current: 234.56608948292705 is not greater than 234.56963579667996. Delta is -0.0035463137529063715 Model is not converging. Current: 327.100573376318 is not greater than 327.10142199921546. Delta is -0.0008486228974788901
Loaded 13 assets, 3404 obs, features=['x_ewvol', 'x_pc1']
Model is not converging. Current: 427.3183204804763 is not greater than 427.32060493028047. Delta is -0.0022844498041649786 Model is not converging. Current: 428.2749629210059 is not greater than 428.2769174246335. Delta is -0.0019545036275872008 Model is not converging. Current: 509.2483583552393 is not greater than 509.24885750700827. Delta is -0.0004991517689632019 Model is not converging. Current: 512.7352541672615 is not greater than 512.7356723270244. Delta is -0.0004181597629440148 Model is not converging. Current: 516.4769012624967 is not greater than 516.4782103939042. Delta is -0.0013091314075381888 Model is not converging. Current: 475.9218608214578 is not greater than 475.9243333213066. Delta is -0.0024724998488068195 Model is not converging. Current: 543.2106435761854 is not greater than 543.2143222719205. Delta is -0.003678695735175097 Model is not converging. Current: 612.2601918751456 is not greater than 612.2668734441668. Delta is -0.006681569021225187 Model is not converging. Current: 606.8076075184647 is not greater than 606.8176610287069. Delta is -0.010053510242187258 Model is not converging. Current: 601.3155875455163 is not greater than 601.321173603596. Delta is -0.005586058079643408 Model is not converging. Current: 595.530505498803 is not greater than 595.5430063848964. Delta is -0.012500886093448571 Model is not converging. Current: 595.4187198374369 is not greater than 595.4430166505438. Delta is -0.02429681310684373 Model is not converging. Current: 780.4332295672039 is not greater than 780.442983818365. Delta is -0.009754251161098182 Model is not converging. Current: 783.4636436217628 is not greater than 783.4762710259783. Delta is -0.012627404215550087 Model is not converging. Current: 661.4717871464429 is not greater than 661.480579846008. Delta is -0.008792699565105977 Model is not converging. Current: 430.485271683611 is not greater than 430.48576104116967. Delta is -0.0004893575586493171 Model is not converging. Current: 414.9587845668064 is not greater than 414.95993299754883. Delta is -0.001148430742432538 Model is not converging. Current: 439.9638217529458 is not greater than 439.96421214238666. Delta is -0.00039038944083813476 Model is not converging. Current: 697.6740478631006 is not greater than 697.7169317616757. Delta is -0.0428838985751554 Model is not converging. Current: 713.3632991821065 is not greater than 713.3677052139906. Delta is -0.004406031884059303 Model is not converging. Current: 702.9921798514354 is not greater than 702.9967879565751. Delta is -0.004608105139709551 Model is not converging. Current: 794.1431466841034 is not greater than 794.1446038038556. Delta is -0.0014571197522172952 Model is not converging. Current: 794.1723698299093 is not greater than 794.1813470361395. Delta is -0.008977206230269985 Model is not converging. Current: 793.4660239450867 is not greater than 793.4720438828839. Delta is -0.006019937797191233 Model is not converging. Current: 794.8513919590821 is not greater than 794.8717084851203. Delta is -0.020316526038186566 Model is not converging. Current: 885.1221199364642 is not greater than 885.1643689885387. Delta is -0.04224905207456686 Model is not converging. Current: 890.1465025362984 is not greater than 890.1794216747718. Delta is -0.03291913847340311 Model is not converging. Current: 893.5227171578038 is not greater than 893.537117674615. Delta is -0.014400516811292619 Model is not converging. Current: 911.049662573319 is not greater than 911.0583039213041. Delta is -0.008641347985076209 Model is not converging. Current: 789.4202157468898 is not greater than 789.4224889140556. Delta is -0.0022731671658675623 Model is not converging. Current: 761.2092350754426 is not greater than 761.2118681294932. Delta is -0.002633054050534156 Model is not converging. Current: 766.7599430281199 is not greater than 766.7630284363755. Delta is -0.003085408255628863
=== Model 3 KPIs (native, pre-cost) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day
Equal-Weight 0.0609 0.0725 0.8405 -0.2009 0.0108
Static HRP 0.0176 0.0198 0.8869 -0.0569 0.0031
IVP (Const-Vol RP, rolling) 0.0172 0.0120 1.4308 -0.0306 0.0017
RA-HRP (prob-blend) 0.0140 0.0097 1.4419 -0.0306 0.0011
RA-HRP (mixture-cov HRP) 0.0148 0.0096 1.5317 -0.0306 0.0010
=== Model 3 KPIs (risk-matched @ 10%) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day Alpha(scale)
Equal-Weight 0.0835 0.1000 0.8351 -0.2683 0.0150 1.3802
Static HRP 0.0875 0.1000 0.8750 -0.2628 0.0155 5.0511
IVP (Const-Vol RP, rolling) 0.1473 0.1000 1.4733 -0.2303 0.0142 8.3250
RA-HRP (prob-blend) 0.1488 0.1000 1.4877 -0.2777 0.0113 10.3164
RA-HRP (mixture-cov HRP) 0.1590 0.1000 1.5898 -0.2793 0.0108 10.3824
Saved net-of-costs grid → ./out/m3_kpi_matchedVol_costGrid.csv
Figures & tables saved in: /Users/jeremy.duriez/Desktop/recherche2/out3
IA + HMM (Neural prior fused with HMM)¶
#!/usr/bin/env python3
# Test 8 — IA + HMM (Neural Prior fused with HMM one-step probabilities) — FAIR VERSION
# --------------------------------------------------------------------------------------
# Core unchanged:
# - K=3, ROLL_FIT_N=500, REFIT_EVERY=3
# - Features x = [log(EW vol, 21d), log(PC1 variance share, 63d)]
# - HMM on x; Tiny MLP learns P_NN(s_t|x_t) on last train window; Fusion: tempered PoE
# - Allocation: probabilistic blend of per-regime MV weights from returns
#
# Fairness/robustness (like Models 1/2/3):
# - Add IVP baseline, ES95 metric, risk-matched (10% vol) headline table
# - Net-of-costs grid (5/10/25/50 bps), subperiod KPIs, turnover stats
# - SAFE turnover density (no SciPy KDE)
# - HMM calibration (Brier/ECE) on fused one-step-ahead probs
# - Data manifest
#
# 10% pb fixed for Static MV:
# - Static MV weights trained once on initial window (no look-ahead), then held
# - When risk-matching Static MV: cap leverage (ALPHA_CAP_STATIC) and use vol floor
# --------------------------------------------------------------------------------------
import os, glob, warnings
warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
DATA_DIR = "./data2"
OUT_DIR = "./out4"
K = 3
ROLL_VOL_N = 21
ROLL_PC_WIN = 63
ROLL_FIT_N = 500
REFIT_EVERY = 3
RIDGE = 1e-6
CLIP = 0.5
SEED = 7
# IA (NN)
FUSE_BETA = 0.40
NN_HIDDEN = 16
NN_EPOCHS = 80
NN_LR = 5e-3
NN_WEIGHT_DECAY = 1e-4
NN_DROPOUT = 0.10
# Fairness/reporting
TARGET_VOL = 0.10
COST_BPS_GRID = [5,10,25,50]
IVP_WIN = 63
IVP_REBAL = 21
# Static baseline safeguards (fix the 10% pb)
STATIC_TRAIN_WIN = 500
ALPHA_CAP_STATIC = 5.0 # <= 5x scaling max when risk-matching static
VOL_FLOOR_STATIC = 0.02 # 2% annual vol floor for static risk-matching
np.random.seed(SEED)
os.makedirs(OUT_DIR, exist_ok=True)
# ------------ Imports (hmmlearn + torch) -------------
try:
from hmmlearn.hmm import GaussianHMM
except Exception as e:
raise SystemExit("pip install hmmlearn\nOriginal error: %r" % e)
try:
import torch
import torch.nn as nn
import torch.optim as optim
except Exception as e:
raise SystemExit("pip install torch\nOriginal error: %r" % e)
DEVICE = "cuda" if torch.cuda.is_available() else "cpu"
torch.manual_seed(SEED)
# ---------------- Helpers ----------------
def load_data_from_dir(path):
files = sorted(glob.glob(os.path.join(path, "*.csv")))
if not files:
raise FileNotFoundError(f"No CSVs found in {path}")
frames = []
for f in files:
tkr = os.path.splitext(os.path.basename(f))[0].upper()
df = pd.read_csv(f)
if "Date" not in df.columns and "date" in df.columns:
df.rename(columns={"date":"Date"}, inplace=True)
close_col = None
for c in ["Close","Adj Close","close","adj_close","AdjClose"]:
if c in df.columns:
close_col = c; break
if close_col is None:
continue
df = df[["Date", close_col]].rename(columns={close_col: tkr})
df["Date"] = pd.to_datetime(df["Date"])
frames.append(df)
if not frames:
raise ValueError("No usable files with Close/Adj Close.")
out = frames[0]
for df in frames[1:]:
out = out.merge(df, on="Date", how="outer")
out.sort_values("Date", inplace=True)
out.set_index("Date", inplace=True)
return out.ffill().dropna(how="all")
def to_log_returns(P): return np.log(P/P.shift(1))
def ew_vol_feature(R, n=21):
ew = R.mean(axis=1)
vol = ew.rolling(n).std()
return np.log(vol + 1e-12).rename("x_ewvol")
def pca_pc1_var_feature(R, win=63):
from numpy.linalg import eigh
vals, idxs = [], []
for i in range(win-1, len(R)):
X = R.iloc[i-win+1:i+1].dropna(how="any")
if len(X) < win:
continue
C = X.cov().values
evals, _ = eigh(C) # ascending
lam1 = float(evals[-1]); total = float(evals.sum()) + 1e-12
vals.append(np.log(lam1/total + 1e-12)); idxs.append(R.index[i])
return pd.Series(vals, index=idxs, name="x_pc1")
def ann_kpis_from_returns(rets, freq=252):
rets = pd.Series(rets).dropna()
if rets.empty: return dict(ret=np.nan, vol=np.nan, sharpe=np.nan, maxdd=np.nan, es95=np.nan)
curve=(1+rets).cumprod()
maxdd=(curve/curve.cummax()-1).min()
ann_ret=(1+rets).prod()**(freq/len(rets)) - 1
ann_vol=rets.std()*np.sqrt(freq)
sharpe=ann_ret/ann_vol if ann_vol>0 else np.nan
q=np.quantile(rets, 0.05)
es95 = -rets[rets<=q].mean() if (rets<=q).any() else np.nan
return dict(ret=ann_ret, vol=ann_vol, sharpe=sharpe, maxdd=maxdd, es95=es95)
def smooth_drawdown_from_returns(rets):
curve=(1+pd.Series(rets).fillna(0)).cumprod()
dd=curve/curve.cummax()-1.0
return curve, dd
def mean_var_weights(mu, Sigma, ridge=1e-6, clip=0.5):
n = len(mu)
S = Sigma + ridge*np.eye(n)
w = np.linalg.pinv(S).dot(mu.reshape(-1,1)).ravel()
w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones(n)/n
w = np.clip(w, -clip, clip)
w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones(n)/n
return w
def inverse_vol_weights(R_window: pd.DataFrame):
vol = R_window.std().replace(0, np.nan)
w = 1.0/vol
w = w.replace([np.inf,-np.inf], np.nan).fillna(0.0)
if w.sum()>0: w = w/w.sum()
else: w = pd.Series(np.ones(len(R_window.columns))/len(R_window.columns), index=R_window.columns)
return w.values
def turnover_series(W: pd.DataFrame):
Wf=W.fillna(method="ffill").fillna(0.0)
dW=(Wf - Wf.shift(1)).abs().sum(axis=1)*0.5
dW.iloc[0]=0.0
return dW
def apply_costs(gross: pd.Series, turnover: pd.Series, cost_bps: float):
c=cost_bps/10000.0
return gross - c*turnover
def risk_match(rets: pd.Series, target_vol=0.10, freq=252, alpha_max=np.inf, vol_floor=0.0):
vol = rets.std()*np.sqrt(freq)
if not np.isfinite(vol) or vol<=0: return rets.copy(), 1.0
denom = max(vol, vol_floor)
alpha = target_vol / denom
if np.isfinite(alpha_max): alpha = min(alpha, alpha_max)
return rets*alpha, alpha
# ---------------- IA (tiny MLP) ----------------
class TinyMLP(nn.Module):
def __init__(self, in_dim, out_dim, hid=16, p=0.1):
super().__init__()
self.net = nn.Sequential(
nn.Linear(in_dim, hid), nn.ReLU(),
nn.Dropout(p),
nn.Linear(hid, hid), nn.ReLU(),
nn.Dropout(p),
nn.Linear(hid, out_dim)
)
def forward(self, x): # logits
return self.net(x)
@torch.no_grad()
def nn_predict_proba(model, X):
model.eval()
logits = model(torch.as_tensor(X, dtype=torch.float32, device=DEVICE))
probs = torch.softmax(logits, dim=-1).cpu().numpy()
return probs
def nn_train(model, X, y, epochs=80, lr=5e-3, wd=1e-4):
model.train()
X_t = torch.as_tensor(X, dtype=torch.float32, device=DEVICE)
y_t = torch.as_tensor(y, dtype=torch.long, device=DEVICE)
opt = optim.Adam(model.parameters(), lr=lr, weight_decay=wd)
loss_fn = nn.CrossEntropyLoss()
for _ in range(epochs):
opt.zero_grad()
logits = model(X_t)
loss = loss_fn(logits, y_t)
loss.backward()
opt.step()
# ---------------- Load data & features ----------------
prices = load_data_from_dir(DATA_DIR)
R = to_log_returns(prices).dropna()
N = R.shape[1]
x_ew = ew_vol_feature(R, n=ROLL_VOL_N)
x_pc = pca_pc1_var_feature(R, win=ROLL_PC_WIN)
x_df = pd.concat([x_ew, x_pc], axis=1).dropna()
R_aligned = R.loc[x_df.index]
dates = x_df.index
tickers = list(R_aligned.columns)
print(f"Loaded {N} assets; feature rows={len(dates)}; features={list(x_df.columns)}")
# Manifest
manifest = pd.DataFrame({
"Ticker": R_aligned.columns,
"FirstDate": [R_aligned.index.min()]*N,
"LastDate": [R_aligned.index.max()]*N,
"Obs": [R_aligned[c].dropna().shape[0] for c in R_aligned.columns]
})
manifest.to_csv(os.path.join(OUT_DIR, "t8_manifest.csv"), index=False)
# ---------------- Rolling IA+HMM ----------------
labels = pd.Series(index=dates, dtype=int)
filt_prob = pd.DataFrame(index=dates, columns=[f"p{k}" for k in range(K)], dtype=float)
ahead_prob = filt_prob.copy()
fused_prob = filt_prob.copy()
model_hmm = None
last_fit = None
# IA model (2→K)
nn_model = TinyMLP(in_dim=x_df.shape[1], out_dim=K, hid=NN_HIDDEN, p=NN_DROPOUT).to(DEVICE)
for i, dt in enumerate(dates):
refit = (last_fit is None) or ((i - last_fit) >= REFIT_EVERY)
# HMM (train / refit on last 500 obs of features)
if refit:
start_fit = max(0, i - ROLL_FIT_N)
X_fit = x_df.iloc[start_fit:i+1].values
if len(X_fit) < max(120, ROLL_FIT_N//4):
continue
model_hmm = GaussianHMM(n_components=K, covariance_type='diag', n_iter=200, random_state=SEED)
model_hmm.fit(X_fit)
last_fit = i
# Obtain teacher labels over the training window for NN (Viterbi path)
path_fit = model_hmm.predict(X_fit)
# Train NN to predict s_t from x_t
try:
nn_train(nn_model, X_fit, path_fit, epochs=NN_EPOCHS, lr=NN_LR, wd=NN_WEIGHT_DECAY)
except Exception as e:
print("NN training skipped due to:", e)
# Inference on rolling window (for probabilities)
start = max(0, i - ROLL_FIT_N)
X_win = x_df.iloc[start:i+1].values
if len(X_win) < 2:
continue
path = model_hmm.predict(X_win)
post = model_hmm.predict_proba(X_win)
labels.iloc[i] = path[-1]
alpha_t = post[-1]
filt_prob.iloc[i] = alpha_t
p_hmm_next = alpha_t @ model_hmm.transmat_
ahead_prob.iloc[i] = p_hmm_next
# IA prior for s_{t+1}: predict s_t from x_t, then push through A
p_nn_t = nn_predict_proba(nn_model, X_win[-1:].copy())[0] # P_NN(s_t | x_t)
p_nn_next = p_nn_t @ model_hmm.transmat_
# Fuse (tempered product of experts)
eps = 1e-12
fused = (p_hmm_next + eps)**(1.0 - FUSE_BETA) * (p_nn_next + eps)**(FUSE_BETA)
fused = fused / fused.sum()
fused_prob.iloc[i] = fused
# ---- Remap states calm→mid→stress by feature stress score
stress_score = (x_df["x_ewvol"] - x_df["x_ewvol"].mean())/(x_df["x_ewvol"].std()+1e-12) + \
(x_df["x_pc1"] - x_df["x_pc1"].mean()) /(x_df["x_pc1"].std() +1e-12)
means = [(k, stress_score[labels==k].mean()) for k in range(K)]
order = [k for (k, _) in sorted(means, key=lambda z:z[1])]
remap = {old:new for new, old in enumerate(order)}
labels = labels.map(remap)
filt_prob = filt_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
ahead_prob = ahead_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
fused_prob = fused_prob.rename(columns={f"p{old}":f"p{remap[old]}" for old in range(K)})
# ---------------- Per-regime weights (MV, NON-CAUSAL, like others) ----------------
def regime_weights_from_labels(R, labels, ridge=1e-6, clip=0.5):
N = R.shape[1]; out = {}
for k in sorted(labels.dropna().unique()):
idx = labels.index[labels==k]
if len(idx) < max(60, N*5):
out[int(k)] = np.ones(N)/N
else:
sub = R.loc[idx]
mu = sub.mean().values
Sigma = sub.cov().values
out[int(k)] = mean_var_weights(mu, Sigma, ridge=ridge, clip=clip)
return out
regime_w = regime_weights_from_labels(R_aligned, labels, RIDGE, CLIP)
# ---------------- Portfolio from fused probabilities ---------------
W_fused = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for dt in R_aligned.index:
p_row = fused_prob.loc[dt]
if p_row.isna().any():
W_fused.loc[dt] = np.ones(N)/N
continue
w = np.zeros(N)
for k in range(K):
w += float(p_row[f"p{k}"]) * regime_w.get(k, np.ones(N)/N)
w = np.clip(w, -CLIP, CLIP)
w = w/np.abs(w).sum() if np.abs(w).sum()>0 else np.ones(N)/N
W_fused.loc[dt] = w
ret_fused = (R_aligned * W_fused.shift(1)).sum(axis=1).fillna(0.0)
turn_fused = turnover_series(W_fused)
# ---------------- Baselines ----------------
# Equal-Weight
w_ew = np.ones(N)/N
W_ew = pd.DataFrame([w_ew]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
ret_ew = (R_aligned @ pd.Series(w_ew, index=R_aligned.columns)).fillna(0.0)
turn_ew = turnover_series(W_ew)
# Static MV — FIXED (train once on initial window; held constant)
train_n = min(STATIC_TRAIN_WIN, len(R_aligned))
Sigma_train = R_aligned.iloc[:train_n].cov().values
mu_train = R_aligned.iloc[:train_n].mean().values
w_static = mean_var_weights(mu_train, Sigma_train, ridge=RIDGE, clip=CLIP)
W_static = pd.DataFrame([w_static]*len(R_aligned), index=R_aligned.index, columns=R_aligned.columns)
ret_static = (R_aligned @ pd.Series(w_static, index=R_aligned.columns)).fillna(0.0)
turn_static = turnover_series(W_static)
# Constant-Vol Risk-Parity (IVP), rolling
W_ivp = pd.DataFrame(index=R_aligned.index, columns=R_aligned.columns, dtype=float)
for i, dt in enumerate(R_aligned.index):
if i<IVP_WIN:
W_ivp.iloc[i]=np.ones(N)/N
else:
if i % IVP_REBAL == 0 or i == IVP_WIN:
Rw=R_aligned.iloc[i-IVP_WIN:i]
W_ivp.iloc[i]=inverse_vol_weights(Rw)
else:
W_ivp.iloc[i]=W_ivp.iloc[i-1]
W_ivp=W_ivp.fillna(method="ffill").fillna(1.0/N)
ret_ivp=(R_aligned * W_ivp.shift(1)).sum(axis=1).fillna(0.0)
turn_ivp = turnover_series(W_ivp)
# ---------------- KPIs (native, pre-cost) ----------------
def krow(name, rets):
k=ann_kpis_from_returns(rets)
return [name, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']]
kpi_native = pd.DataFrame(
[krow("Equal-Weight", ret_ew),
krow("Static MV", ret_static),
krow("IVP (Const-Vol RP, rolling)", ret_ivp),
krow("IA+HMM (fused)", ret_fused)],
columns=["Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"]
)
kpi_native.to_csv(os.path.join(OUT_DIR,"t8_kpi_native.csv"), index=False)
print("\n=== Test 8 KPIs (native, pre-cost) ===")
print(kpi_native.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
# ---------------- Risk-matched @ 10% (with static cap/floor) ---------------
def matched_series_and_alpha(name, rets, turn, *, alpha_max=np.inf, vol_floor=0.0):
rm, alpha = risk_match(rets, TARGET_VOL, alpha_max=alpha_max, vol_floor=vol_floor)
return name, rm, turn*alpha, alpha, ann_kpis_from_returns(rm)
matched = {}
turn_m = {}
alphas = {}
rows=[]
for name, rets, turn in [
("Equal-Weight", ret_ew, turn_ew),
("Static MV", ret_static, turn_static),
("IVP (Const-Vol RP, rolling)", ret_ivp, turn_ivp),
("IA+HMM (fused)", ret_fused, turn_fused),
]:
if name == "Static MV":
nm, rm, tm, a, k = matched_series_and_alpha(name, rets, turn,
alpha_max=ALPHA_CAP_STATIC,
vol_floor=VOL_FLOOR_STATIC)
else:
nm, rm, tm, a, k = matched_series_and_alpha(name, rets, turn)
matched[nm]=rm; turn_m[nm]=tm; alphas[nm]=a
rows.append([nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95'], a])
kpi_matched = pd.DataFrame(rows, columns=["Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day","Alpha(scale)"])
kpi_matched.to_csv(os.path.join(OUT_DIR,"t8_kpi_matchedVol.csv"), index=False)
print("\n=== Test 8 KPIs (risk-matched @ {:.0f}%) ===".format(TARGET_VOL*100))
print(kpi_matched.to_string(index=False, float_format=lambda x: f"{x:,.4f}"))
# ---------------- Net-of-costs grid ---------------------------------------
def cost_grid(cost_list):
out=[]
for c in cost_list:
for nm in matched:
net = apply_costs(matched[nm], turn_m[nm], c)
k = ann_kpis_from_returns(net)
out.append([c, nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
return pd.DataFrame(out, columns=["Cost_bps","Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"])
kpi_costs = cost_grid(COST_BPS_GRID)
kpi_costs.to_csv(os.path.join(OUT_DIR,"t8_kpi_matchedVol_costGrid.csv"), index=False)
print("\nSaved net-of-costs grid → ./out/t8_kpi_matchedVol_costGrid.csv")
# ---------------- Turnover stats & SAFE density plot -----------------------
turn_stats=[]
for nm, ts in turn_m.items():
turn_stats.append([nm, ts.mean(), ts.median(), ts.quantile(0.9), ts.max()])
turn_df = pd.DataFrame(turn_stats, columns=["Strategy","Turnover_mean","Turnover_median","Turnover_p90","Turnover_max"])
turn_df.to_csv(os.path.join(OUT_DIR,"t8_turnover_stats.csv"), index=False)
def plot_turnover_distribution(turn_map, out_path):
fig, ax = plt.subplots(figsize=(10,4))
for name, ts in turn_map.items():
s = pd.Series(ts).replace([np.inf, -np.inf], np.nan).dropna()
s = s[np.isfinite(s)]
if (len(s) == 0) or (s.std() <= 1e-12) or (s.nunique() < 3):
val = float(s.iloc[0]) if len(s) else 0.0
ax.axvline(val, linestyle="--", linewidth=1.5, label=f"{name} (spike)")
else:
counts, bins = np.histogram(s, bins=40, density=True)
centers = 0.5*(bins[1:] + bins[:-1])
ax.plot(centers, counts, linewidth=1.5, label=name)
ax.set_title("Turnover distribution (one-way, risk-matched)")
ax.grid(True); ax.legend()
fig.tight_layout(); fig.savefig(out_path, dpi=180)
plot_turnover_distribution(turn_m, os.path.join(OUT_DIR,"t8_fig_turnover_density.png"))
# ---------------- Subperiod KPIs (risk-matched) ----------------------------
def subperiod_slices(idx):
eras=[("2000-01-01","2009-12-31"), ("2010-01-01","2019-12-31"), ("2020-01-01", str(idx.max().date()))]
out=[]
for s,e in eras:
sdt=pd.Timestamp(s); edt=pd.Timestamp(e)
mask=(idx>=sdt)&(idx<=edt)
if mask.any(): out.append((s,e,mask))
return out
subs=[]
for (s,e,mask) in subperiod_slices(R_aligned.index):
for nm in matched:
k=ann_kpis_from_returns(matched[nm].loc[mask])
subs.append([f"{s}–{e}", nm, k['ret'], k['vol'], k['sharpe'], k['maxdd'], k['es95']])
sub_df=pd.DataFrame(subs, columns=["Period","Strategy","Ann.Return","Ann.Vol","Sharpe","MaxDD","ES95_day"])
sub_df.to_csv(os.path.join(OUT_DIR,"t8_kpi_matchedVol_subperiods.csv"), index=False)
# ---------------- HMM calibration (one-step-ahead, fused) ------------------
probs = fused_prob.copy().dropna()
y_true = labels.reindex(probs.index).shift(-1).dropna()
probs = probs.loc[y_true.index]
Y = np.zeros((len(y_true), K))
for i, y in enumerate(y_true.values.astype(int)):
if 0 <= y < K: Y[i, y]=1.0
P = probs[[f"p{k}" for k in range(K)]].values
brier = np.mean(np.sum((P - Y)**2, axis=1))
conf = P.max(axis=1); pred = P.argmax(axis=1); acc = (pred==y_true.values.astype(int)).astype(float)
bins = np.linspace(0,1,11); ece=0.0
for b0, b1 in zip(bins[:-1], bins[1:]):
idx = (conf>=b0) & (conf<b1)
if idx.any():
ece += np.mean(idx) * abs(acc[idx].mean() - conf[idx].mean())
pd.DataFrame({"Metric":["Brier","ECE_top10"], "Value":[brier, ece]}).to_csv(os.path.join(OUT_DIR,"t8_hmm_calibration.csv"), index=False)
# ---------------- Figures ---------------------------------------------------
# Risk-matched log cumulative wealth (pre-cost)
plt.figure(figsize=(10,4))
for nm, sr in matched.items():
curve=(1+sr).cumprod()
plt.plot(curve.index, np.log(curve), label=nm)
plt.title("Log Cumulative Wealth — Risk-matched (10% vol)") # Static may be <10% if capped
plt.legend(); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"t8_fig_cum_wealth_matched.png"), dpi=180)
# Risk-matched drawdowns
def ddplot_from_rets(name, rets):
curve, dd = smooth_drawdown_from_returns(rets)
plt.figure(figsize=(10,3)); plt.plot(dd)
plt.title(f"Drawdown (risk-matched) — {name}"); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR, f"t8_fig_drawdown_matched_{name.replace(' ','_').replace('-','')}.png"), dpi=180)
for nm, sr in matched.items():
ddplot_from_rets(nm, sr)
# Native wealth (pre-cost)
curve_fused=(1+ret_fused).cumprod(); curve_static=(1+ret_static).cumprod()
curve_ivp=(1+ret_ivp).cumprod(); curve_ew=(1+ret_ew).cumprod()
plt.figure(figsize=(10,4))
plt.plot(np.log(curve_fused), label="IA+HMM (fused)")
plt.plot(np.log(curve_static), label="Static MV")
plt.plot(np.log(curve_ivp), label="IVP (rolling)")
plt.plot(np.log(curve_ew), label="Equal-Weight")
plt.title("Log Cumulative Wealth — Test 8 (native)"); plt.legend(); plt.grid(True, alpha=.3); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"t8_cum_wealth_native.png"), dpi=180)
# Probability ribbon (fused)
cols=[f"p{k}" for k in range(K)]
pdf=fused_prob[cols].dropna(how="any")
plt.figure(figsize=(12,3.6))
plt.stackplot(pdf.index, pdf.values.T, labels=cols, alpha=0.9)
plt.ylim(0,1); plt.title("One-step-ahead Probabilities — IA+HMM (fused)"); plt.legend(ncol=K, fontsize=8, frameon=False)
plt.grid(True, alpha=.2); plt.tight_layout()
plt.savefig(os.path.join(OUT_DIR,"t8_prob_ribbon.png"), dpi=180)
# Weights timeline (fused)
plt.figure(figsize=(11,4))
plt.stackplot(W_fused.index, W_fused.fillna(method="ffill").values.T, labels=tickers)
plt.title("Portfolio Weights — IA+HMM (fused)"); plt.legend(loc="upper left", ncol=min(N,4), fontsize=8)
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR,"t8_weights_fused.png"), dpi=180)
print(f"\nFigures & tables saved in: {os.path.abspath(OUT_DIR)}")
Model is not converging. Current: 389.87283970835267 is not greater than 389.91854359257286. Delta is -0.04570388422018823
Loaded 13 assets; feature rows=3362; features=['x_ewvol', 'x_pc1']
Model is not converging. Current: 395.60349779898564 is not greater than 395.6122636760006. Delta is -0.008765877014980106 Model is not converging. Current: 404.2599203217861 is not greater than 404.26936548010417. Delta is -0.009445158318044378 Model is not converging. Current: 420.36977301248635 is not greater than 420.38081938787917. Delta is -0.011046375392822938 Model is not converging. Current: 532.7838403663428 is not greater than 532.8045589952701. Delta is -0.020718628927284044 Model is not converging. Current: 523.2122200256835 is not greater than 523.2601395497974. Delta is -0.04791952411392231 Model is not converging. Current: 709.8789261884571 is not greater than 709.8798558887838. Delta is -0.0009297003267647597 Model is not converging. Current: 692.1376087245156 is not greater than 692.1396690939667. Delta is -0.0020603694511009962 Model is not converging. Current: 677.1086596972053 is not greater than 677.1110946649715. Delta is -0.0024349677661348323 Model is not converging. Current: 497.3941108311458 is not greater than 497.3941870640186. Delta is -7.623287280011937e-05 Model is not converging. Current: 647.2603760262214 is not greater than 647.2612236604854. Delta is -0.0008476342640051371 Model is not converging. Current: 638.4778442645846 is not greater than 638.4848189671235. Delta is -0.0069747025388551265 Model is not converging. Current: 497.62463232852065 is not greater than 497.62496636703685. Delta is -0.00033403851620050773 Model is not converging. Current: 468.37887909742966 is not greater than 468.3835790480223. Delta is -0.004699950592623736 Model is not converging. Current: 473.0591907627078 is not greater than 473.0661737642574. Delta is -0.006983001549656365 Model is not converging. Current: 560.0942112781146 is not greater than 560.0944139239596. Delta is -0.000202645844979088 Model is not converging. Current: 561.9855916122344 is not greater than 561.9859251994662. Delta is -0.0003335872318075417 Model is not converging. Current: 560.9462865959019 is not greater than 560.9467565087089. Delta is -0.00046991280692054715 Model is not converging. Current: 569.9735594618257 is not greater than 569.9888135330621. Delta is -0.015254071236313393 Model is not converging. Current: 513.9888519080189 is not greater than 514.017249539479. Delta is -0.02839763146016594 Model is not converging. Current: 505.9065507840104 is not greater than 505.9325281762053. Delta is -0.02597739219493178 Model is not converging. Current: 496.8747461679629 is not greater than 496.8940313340508. Delta is -0.01928516608791142 Model is not converging. Current: 493.10471063548647 is not greater than 493.1100064974062. Delta is -0.005295861919705658 Model is not converging. Current: 453.5731327851453 is not greater than 453.5732345809892. Delta is -0.00010179584393199548 Model is not converging. Current: 444.7933227674599 is not greater than 444.7935453070773. Delta is -0.00022253961736851124 Model is not converging. Current: 398.01333655053435 is not greater than 398.0281765801007. Delta is -0.01484002956635777 Model is not converging. Current: 390.7892675515517 is not greater than 390.80067757864106. Delta is -0.011410027089368668 Model is not converging. Current: 406.3594393384705 is not greater than 406.38140047739375. Delta is -0.02196113892324547 Model is not converging. Current: 408.57293892512826 is not greater than 408.59163249167705. Delta is -0.01869356654879084 Model is not converging. Current: 422.517630334027 is not greater than 422.5401400854047. Delta is -0.022509751377697285 Model is not converging. Current: 427.84845482848306 is not greater than 427.87875449935547. Delta is -0.03029967087240948 Model is not converging. Current: 432.9376840215134 is not greater than 432.93894313415836. Delta is -0.0012591126449592593 Model is not converging. Current: 429.7194932000672 is not greater than 429.7201140128609. Delta is -0.0006208127937270547 Model is not converging. Current: 444.63636410513374 is not greater than 444.67077087383865. Delta is -0.03440676870491188 Model is not converging. Current: 458.890811536571 is not greater than 458.9368395188587. Delta is -0.04602798228768279 Model is not converging. Current: 501.7909691715577 is not greater than 501.80617906783107. Delta is -0.015209896273347567 Model is not converging. Current: 516.968558607946 is not greater than 516.9793435893532. Delta is -0.010784981407255145 Model is not converging. Current: 566.9588058217353 is not greater than 566.959076575291. Delta is -0.0002707535556965013 Model is not converging. Current: 576.9378295435155 is not greater than 576.9395475683427. Delta is -0.001718024827255249 Model is not converging. Current: 563.3181894420273 is not greater than 563.3188705877014. Delta is -0.0006811456740933863 Model is not converging. Current: 566.5981866707649 is not greater than 566.5982157718776. Delta is -2.910111265919113e-05 Model is not converging. Current: 559.3862085574262 is not greater than 559.3910065990974. Delta is -0.004798041671165265 Model is not converging. Current: 565.6964827430708 is not greater than 565.6999572398837. Delta is -0.003474496812941652 Model is not converging. Current: 562.0576815375418 is not greater than 562.0609985615853. Delta is -0.003317024043440142 Model is not converging. Current: 571.18586272219 is not greater than 571.1881566951915. Delta is -0.0022939730015423265 Model is not converging. Current: 559.1272190083262 is not greater than 559.1291104979588. Delta is -0.0018914896326123198 Model is not converging. Current: 602.2516234704003 is not greater than 602.2524558736724. Delta is -0.000832403272056581 Model is not converging. Current: 623.3517198848233 is not greater than 623.3529717342053. Delta is -0.0012518493820152798 Model is not converging. Current: 643.2416420327299 is not greater than 643.2429544022963. Delta is -0.001312369566335292 Model is not converging. Current: 645.9588758596371 is not greater than 645.9594446361631. Delta is -0.0005687765259381194 Model is not converging. Current: 641.725184585292 is not greater than 641.7292498974475. Delta is -0.004065312155489664 Model is not converging. Current: 640.6099708196236 is not greater than 640.6100062032195. Delta is -3.5383595900384535e-05 Model is not converging. Current: 641.1182704187269 is not greater than 641.121223220791. Delta is -0.002952802064100979 Model is not converging. Current: 642.4738808363991 is not greater than 642.4779252570721. Delta is -0.004044420673039895 Model is not converging. Current: 646.2972767632926 is not greater than 646.3010637404075. Delta is -0.003786977114828005 Model is not converging. Current: 648.2628522276609 is not greater than 648.2657217857371. Delta is -0.002869558076213252 Model is not converging. Current: 781.17901881952 is not greater than 781.1807263215957. Delta is -0.0017075020756465165 Model is not converging. Current: 844.0410423180136 is not greater than 844.0419943069044. Delta is -0.0009519888907334462 Model is not converging. Current: 821.6122452382398 is not greater than 821.6142392136907. Delta is -0.0019939754508868646 Model is not converging. Current: 608.4144776748001 is not greater than 608.4171884011388. Delta is -0.0027107263387051717 Model is not converging. Current: 615.5738343544865 is not greater than 615.5767204248887. Delta is -0.002886070402155383 Model is not converging. Current: 691.0862571102078 is not greater than 691.0904891170456. Delta is -0.0042320068378103315 Model is not converging. Current: 691.884083607118 is not greater than 691.8883400837645. Delta is -0.004256476646560259 Model is not converging. Current: 691.5933614977193 is not greater than 691.5981597205036. Delta is -0.004798222784302197 Model is not converging. Current: 636.1451718217206 is not greater than 636.1457524068984. Delta is -0.0005805851777722637 Model is not converging. Current: 646.7160351641934 is not greater than 647.1561691657477. Delta is -0.44013400155427007 Model is not converging. Current: 661.8167191709838 is not greater than 662.0372497037188. Delta is -0.22053053273498335 Model is not converging. Current: 667.1582536874924 is not greater than 667.4221345743808. Delta is -0.26388088688838707 Model is not converging. Current: 672.7331673794818 is not greater than 673.0540673213847. Delta is -0.32089994190289417 Model is not converging. Current: 666.7998120075691 is not greater than 667.0102819612648. Delta is -0.21046995369567867 Model is not converging. Current: 666.3059014283275 is not greater than 666.4048955082361. Delta is -0.09899407990860709 Model is not converging. Current: 613.440410328221 is not greater than 613.5257559376993. Delta is -0.0853456094782814 Model is not converging. Current: 606.1417442396211 is not greater than 606.1951246358524. Delta is -0.05338039623131863 Model is not converging. Current: 600.7891245374491 is not greater than 600.8618575142059. Delta is -0.07273297675681079 Model is not converging. Current: 590.0162507323257 is not greater than 590.1353875098237. Delta is -0.1191367774979426 Model is not converging. Current: 606.3707454253064 is not greater than 606.4580413214413. Delta is -0.08729589613494682 Model is not converging. Current: 605.4732142568412 is not greater than 605.5274956160399. Delta is -0.054281359198739665 Model is not converging. Current: 605.5176509042411 is not greater than 605.525586811296. Delta is -0.007935907054843483 Model is not converging. Current: 606.2565081533663 is not greater than 606.2848715353983. Delta is -0.02836338203201194 Model is not converging. Current: 606.7509265233965 is not greater than 606.7553221946939. Delta is -0.004395671297402259 Model is not converging. Current: 558.8773167497272 is not greater than 558.8813221806857. Delta is -0.004005430958500256 Model is not converging. Current: 597.0666982370659 is not greater than 597.0797391331097. Delta is -0.013040896043776229 Model is not converging. Current: 554.752573775024 is not greater than 554.7565105457054. Delta is -0.003936770681434609 Model is not converging. Current: 553.1359039678034 is not greater than 553.1395933998209. Delta is -0.0036894320174951645 Model is not converging. Current: 551.4789157504738 is not greater than 551.5515755562857. Delta is -0.07265980581189524 Model is not converging. Current: 663.3510832251052 is not greater than 663.3841216375271. Delta is -0.0330384124218881 Model is not converging. Current: 656.1958205050584 is not greater than 656.21101485804. Delta is -0.01519435298166627 Model is not converging. Current: 620.2789008812879 is not greater than 620.2820550931145. Delta is -0.0031542118266543184 Model is not converging. Current: 622.9654494902916 is not greater than 622.9906149110335. Delta is -0.025165420741927846 Model is not converging. Current: 664.317944780949 is not greater than 664.3181075029149. Delta is -0.0001627219659212642 Model is not converging. Current: 687.4575429999189 is not greater than 687.4623321099035. Delta is -0.004789109984585593 Model is not converging. Current: 693.5398754709931 is not greater than 693.5629393372019. Delta is -0.023063866208758554 Model is not converging. Current: 685.5607089768239 is not greater than 685.6270437672385. Delta is -0.06633479041465762 Model is not converging. Current: 685.6979303426833 is not greater than 685.7729752831822. Delta is -0.0750449404988558 Model is not converging. Current: 685.307338792074 is not greater than 685.3878031435186. Delta is -0.08046435144456154 Model is not converging. Current: 684.7568930217823 is not greater than 684.8396052391329. Delta is -0.08271221735060408 Model is not converging. Current: 691.64895167386 is not greater than 691.6602250415447. Delta is -0.011273367684680125 Model is not converging. Current: 686.175677832811 is not greater than 686.2529891027104. Delta is -0.07731126989949644 Model is not converging. Current: 683.7305017671897 is not greater than 683.7956175277266. Delta is -0.0651157605368553 Model is not converging. Current: 657.1410721571496 is not greater than 657.1528432752615. Delta is -0.01177111811193754 Model is not converging. Current: 626.5788044697925 is not greater than 626.5833217805884. Delta is -0.00451731079590445 Model is not converging. Current: 609.9742723842236 is not greater than 609.9932929883287. Delta is -0.019020604105094208 Model is not converging. Current: 610.8768527483965 is not greater than 610.8803358151295. Delta is -0.003483066732997031 Model is not converging. Current: 607.7430210659807 is not greater than 607.758386675833. Delta is -0.015365609852324269 Model is not converging. Current: 605.2001156266639 is not greater than 605.2177796053794. Delta is -0.01766397871551817 Model is not converging. Current: 632.2520822580158 is not greater than 632.3096044689308. Delta is -0.05752221091506726 Model is not converging. Current: 624.7246277998114 is not greater than 624.8197662184324. Delta is -0.0951384186209907 Model is not converging. Current: 633.8707062631968 is not greater than 633.8725145960686. Delta is -0.0018083328718603298 Model is not converging. Current: 691.1557900097373 is not greater than 691.1605245976957. Delta is -0.004734587958409975 Model is not converging. Current: 696.3569925712039 is not greater than 696.3576699852701. Delta is -0.0006774140662173522 Model is not converging. Current: 735.9966012205355 is not greater than 736.1497463670136. Delta is -0.15314514647809574 Model is not converging. Current: 741.4090534127449 is not greater than 741.5674348477294. Delta is -0.1583814349844488 Model is not converging. Current: 806.2690004214189 is not greater than 806.2716947495585. Delta is -0.0026943281395688246 Model is not converging. Current: 870.4417880625542 is not greater than 870.4421267040966. Delta is -0.0003386415423847211 Model is not converging. Current: 874.2910441909661 is not greater than 874.2911536886778. Delta is -0.00010949771171908651 Model is not converging. Current: 902.5394478540348 is not greater than 902.5482488208984. Delta is -0.008800966863532267 Model is not converging. Current: 907.9117738031123 is not greater than 907.9208997807042. Delta is -0.00912597759190703 Model is not converging. Current: 917.1479211180495 is not greater than 917.1547300293647. Delta is -0.006808911315260957 Model is not converging. Current: 920.6646704957557 is not greater than 920.6710136945203. Delta is -0.006343198764625413 Model is not converging. Current: 921.6073421545094 is not greater than 921.6138747908119. Delta is -0.006532636302495121 Model is not converging. Current: 923.3575380980798 is not greater than 923.3601991978848. Delta is -0.002661099804981859 Model is not converging. Current: 944.8415612281647 is not greater than 944.8418222631809. Delta is -0.00026103501625129866 Model is not converging. Current: 928.4859346726997 is not greater than 928.4984065529873. Delta is -0.012471880287534987 Model is not converging. Current: 927.6347572647405 is not greater than 927.6385253390098. Delta is -0.003768074269260069 Model is not converging. Current: 935.9690543917699 is not greater than 935.9699588287401. Delta is -0.0009044369702451149 Model is not converging. Current: 944.7283235016441 is not greater than 944.7392140943513. Delta is -0.010890592707141877 Model is not converging. Current: 953.3599412670986 is not greater than 953.4700051737983. Delta is -0.11006390669967914 Model is not converging. Current: 967.4285177642461 is not greater than 967.5282016215652. Delta is -0.09968385731917806 Model is not converging. Current: 970.1763393033673 is not greater than 970.212708906224. Delta is -0.03636960285666646 Model is not converging. Current: 954.5609003640387 is not greater than 954.5762753662549. Delta is -0.015375002216160283 Model is not converging. Current: 1025.9247075782428 is not greater than 1026.023750816467. Delta is -0.0990432382243398 Model is not converging. Current: 1055.636349711419 is not greater than 1055.73155148508. Delta is -0.09520177366107418 Model is not converging. Current: 1047.1504434056817 is not greater than 1047.1900408399013. Delta is -0.039597434219558636 Model is not converging. Current: 1068.7605154622584 is not greater than 1068.901361767094. Delta is -0.14084630483557703 Model is not converging. Current: 1070.4772474350045 is not greater than 1070.6538082224338. Delta is -0.17656078742925274 Model is not converging. Current: 1063.5499616253587 is not greater than 1063.6107377776837. Delta is -0.06077615232493372 Model is not converging. Current: 1100.2480118703263 is not greater than 1100.2560587025926. Delta is -0.008046832266245474 Model is not converging. Current: 1143.2728553065215 is not greater than 1143.3715056965918. Delta is -0.09865039007036103 Model is not converging. Current: 1141.1088328479393 is not greater than 1141.199624631882. Delta is -0.09079178394267728 Model is not converging. Current: 1136.4830069710097 is not greater than 1136.5606864059755. Delta is -0.0776794349658303 Model is not converging. Current: 1132.6185676430289 is not greater than 1132.682428731013. Delta is -0.0638610879841508 Model is not converging. Current: 1128.5292534945431 is not greater than 1128.5758735179184. Delta is -0.04662002337522608 Model is not converging. Current: 1127.056686845815 is not greater than 1127.1168752200665. Delta is -0.060188374251538335 Model is not converging. Current: 1122.2897396306244 is not greater than 1122.3147846895279. Delta is -0.025045058903515383 Model is not converging. Current: 1142.1820087454985 is not greater than 1142.1862650192131. Delta is -0.004256273714645431 Model is not converging. Current: 1099.0578586872014 is not greater than 1099.0620589242747. Delta is -0.004200237073291646 Model is not converging. Current: 1135.6240678363877 is not greater than 1135.6240975103137. Delta is -2.9673926064788247e-05 Model is not converging. Current: 1064.4736230892174 is not greater than 1064.4879871150395. Delta is -0.014364025822033 Model is not converging. Current: 1005.2489913743328 is not greater than 1005.2498028132604. Delta is -0.0008114389275988287 Model is not converging. Current: 1042.2737712841042 is not greater than 1042.2785253362642. Delta is -0.004754052160024003 Model is not converging. Current: 1037.0706277696772 is not greater than 1037.0768267434203. Delta is -0.00619897374303946 Model is not converging. Current: 1034.8494148701352 is not greater than 1034.853882810512. Delta is -0.004467940376798651 Model is not converging. Current: 816.1427352392408 is not greater than 816.1470397044515. Delta is -0.004304465210680064 Model is not converging. Current: 736.8435962831383 is not greater than 736.867376018414. Delta is -0.0237797352757525 Model is not converging. Current: 736.2549722902698 is not greater than 736.2701635859863. Delta is -0.015191295716476816 Model is not converging. Current: 733.0088374746977 is not greater than 733.0445230081223. Delta is -0.035685533424612004 Model is not converging. Current: 730.9667184427168 is not greater than 730.9965400264886. Delta is -0.02982158377176347 Model is not converging. Current: 732.2452023818581 is not greater than 732.2527216315566. Delta is -0.007519249698475505 Model is not converging. Current: 733.5755413108765 is not greater than 733.5862220206391. Delta is -0.010680709762596052 Model is not converging. Current: 734.2816026619587 is not greater than 734.2929153267506. Delta is -0.011312664791944371 Model is not converging. Current: 732.6312220118732 is not greater than 732.6621145275719. Delta is -0.03089251569872431 Model is not converging. Current: 717.5379288239386 is not greater than 717.5592496664099. Delta is -0.021320842471254764 Model is not converging. Current: 710.2769744572315 is not greater than 710.3001240347056. Delta is -0.023149577474100624 Model is not converging. Current: 704.2489738721843 is not greater than 704.2684253688623. Delta is -0.019451496678016156 Model is not converging. Current: 699.0026015397397 is not greater than 699.0122454075054. Delta is -0.009643867765703362 Model is not converging. Current: 694.0550121024608 is not greater than 694.0671688671828. Delta is -0.012156764722021762 Model is not converging. Current: 690.2343320702761 is not greater than 690.2417032675322. Delta is -0.00737119725613411 Model is not converging. Current: 728.8780970422679 is not greater than 728.8802112676173. Delta is -0.0021142253493735552 Model is not converging. Current: 726.6178356581345 is not greater than 726.6183299051031. Delta is -0.0004942469686284312 Model is not converging. Current: 724.6625458449145 is not greater than 724.6626724694527. Delta is -0.00012662453821121744 Model is not converging. Current: 692.7143479235503 is not greater than 692.7326337464873. Delta is -0.018285822936945806
=== Test 8 KPIs (native, pre-cost) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day
Equal-Weight 0.0602 0.0728 0.8267 -0.2053 0.0110
Static MV 0.0197 0.0334 0.5894 -0.0964 0.0050
IVP (Const-Vol RP, rolling) 0.0199 0.0118 1.6913 -0.0216 0.0016
IA+HMM (fused) 0.0282 0.0136 2.0786 -0.0271 0.0017
=== Test 8 KPIs (risk-matched @ 10%) ===
Strategy Ann.Return Ann.Vol Sharpe MaxDD ES95_day Alpha(scale)
Equal-Weight 0.0821 0.1000 0.8211 -0.2728 0.0151 1.3737
Static MV 0.0566 0.1000 0.5657 -0.2683 0.0151 2.9969
IVP (Const-Vol RP, rolling) 0.1771 0.1000 1.7707 -0.1717 0.0136 8.4931
IA+HMM (fused) 0.2222 0.1000 2.2215 -0.1869 0.0129 7.3772
Saved net-of-costs grid → ./out/t8_kpi_matchedVol_costGrid.csv
Figures & tables saved in: /Users/jeremy.duriez/Desktop/recherche2/out4
